System And Method Of Conducting In Silico Clinical Trials

ABSTRACT

Described herein are computer systems and methods for generating, in silico, one or more virtual patients. The one or more virtual patients can be used, for example, to conduct a virtual clinical trial, such as to determine an efficacy of a therapy or medical device. The one or more virtual patients mathematically represent a physiological system in an actual patient for a health condition. Also, the one or more virtual patients can be used, for example, to determine an adverse effect from a therapy or any other deviation from a therapy, or compliance of a patient suffering from a health condition on a therapeutic protocol.

RELATED APPLICATION

The present application claims priority to provisional application No.62/093,977, entitled “System And Method Of Conducting In Silico ClinicalTrials,” filed on Dec. 18, 2014. This application incorporates thecontent of the provisional application by reference in its entirety.

BACKGROUND OF THE INVENTION

Clinical trials are conducted to collect data and other informationregarding the safety and efficacy of drugs and medical devices. Thereare several steps and stages of approval in the clinical trials processbefore a drug or device can be sold in the consumer market. These stepstypically require millions of dollars of investment and take years tocomplete.

Drug and device testing begins with extensive laboratory research whichcan involve years of experiments in animals and human cells. If theinitial laboratory research is successful, the Food and DrugAdministration (FDA) can approve to continue research and testing inhumans. Once approved, human testing of drugs and medical devices canbegin and is typically conducted in four phases. Each phase isconsidered a separate trial and, after completion of a phase,investigators are required to submit their data for approval from theFDA before continuing to the next phase.

A Phase I clinical trial assesses the safety of a drug or device inhealthy volunteers. This initial phase of testing, which can takeseveral months to complete, usually includes a small number of healthyvolunteers (20 to 100). The study is designed to determine the effectsof the drug or device on humans including how it is absorbed,metabolized, and excreted. This phase also investigates the side effectsthat occur as dosage levels are increased.

A Phase II clinical trial tests the efficacy of a drug or medicaldevice. This second phase of testing can last from several months to acouple of years, and typically involves up to several hundred patients.Most phase II studies are randomized trials where one group of patientsreceives the new drug, while a second “control” group receives astandard treatment or placebo. Often these studies are “blinded,”meaning that neither the patients nor the researchers know who hasreceived the drug. This allows investigators to provide thepharmaceutical company and the FDA with comparative information aboutthe relative safety and effectiveness of the new drug.

A Phase III clinical trial involves randomized and blind testing inseveral hundred to several thousand patients. This large-scale testing,which can last several years, provides a more thorough understanding ofthe effectiveness of the drug or device as well as the range of possibleadverse reactions. Once Phase III is complete, a pharmaceutical companycan request FDA approval for marketing the drug.

Finally, a Phase IV clinical trial, often called Post MarketingSurveillance Trial, can be conducted after a drug or device has beenapproved for consumer sale. Objectives at this stage are: (1) to comparea drug with other drugs already in the market; (2) to monitor a drug'slong term effectiveness and impact on a patient's quality of life; and(3) to determine the cost-effectiveness of a drug therapy relative toother traditional and new therapies. Phase IV studies can result in adrug or device being taken off the market or restrictions of use beingplaced on the product depending on the findings in the study.

Conducting all phases of a clinical trial to test a new drug or medicaldevice presents many challenges to investigators. Some challengesinclude the time and financial demands of clinical trials, thecomplexity of FDA regulations, inadequate research training among manyclinicians, pressures to which both investigators and volunteers aresubjected during the course of a trial, and data collection challenges(e.g., medical records, quality control).

Therefore, there is a need for reliable, cost-effective predictivemodels for performing a clinical trial, in silico, to supplement and/orreplace actual clinical trials.

SUMMARY OF THE INVENTION

Described herein are methods and computer systems for generating, insilico, one or more virtual patients. The one or more virtual patientscan be used, for example, to conduct a virtual clinical trial, e.g., todetermine safety and efficacy of a therapy. The virtual patients canmathematically represent one or more physiological systems in an actualpatient. As discussed in more detail below, the virtual patients can beused, for example, to assess the effect of a therapy, to optimize atherapy for administration to actual patients, and to reject a therapybased on observed adverse effects on virtual patients, etc.

In one aspect, the present invention relates to a method for determiningefficacy of a therapy. The method comprises generating in silico aplurality of virtual patients for modeling a health condition based ondata collected from a population of previously treated patients. Thecollected data represents at least one measured biological response ofthe previously treated patients to a previously administered therapeuticregimen. The method further comprises applying a simulated therapy tothe virtual patients and determining one or more physiologicalparameters in the virtual patients in response to the simulated therapy.

The above method for determining efficacy of a therapy can furthercomprise adjusting the simulated therapy based on the determined one ormore physiological parameters of the virtual patients, applying theadjusted simulated therapy to the virtual patients, and determining oneor more physiological parameters in response to the adjusted simulatedtherapy.

In some embodiments, the method further comprises repeating theadjusting, applying, and determining steps so as to optimize thesimulated therapy for application to one or more actual patients.

In a related aspect, the invention relates to a method for performing aclinical trial in silico. The method comprises generating in silico oneor more virtual patients, wherein the one or more virtual patients aregenerated based on data collected from a population of previouslytreated patients suffering from a health condition. The method furthercomprises applying a simulated therapy to each of the one or morevirtual patients and determining one or more physiological parameters inthe virtual patients in response to the simulated therapy.

In some embodiments, the simulated therapy includes a simulatednon-compliance with a prescribed therapy.

In one embodiment, the above method for performing a clinical trial canfurther comprise the following iterative steps to derive an optimaltherapy for the health condition: modifying one or more parameters ofthe simulated therapy, thereby creating a modified simulated therapy,administering the modified simulated therapy to each of the virtualpatients, and determining one or more physiological parameters in thevirtual patients in response to the modified simulated therapy.

In one aspect, the invention relates to a computer system fordetermining efficacy of a therapy. The computer system comprises ageneration module to generate in silico a plurality of virtual patientsfor modeling a health condition based on data collected from apopulation of previously treated patients suffering from the healthcondition, where the collected data represents at least one measuredbiological response of the previously treated patients to a previouslyadministered therapeutic regimen, and a simulation module to apply asimulated therapy to the virtual patients and to determine one or morephysiological parameters of the virtual patients in response to thesimulated therapy.

In one embodiment of the computer system, the simulation module isfurther configured to adjust the simulated therapy based on the one ormore physiological parameters thereby creating a modified simulatedtherapy. The simulation module applies the modified simulated therapy tothe virtual patients, and determines one or more physiologicalparameters of the virtual patients in response to the applied modifiedsimulated therapy.

In one aspect, the invention relates to a method of generating a virtualpatient. The method comprises generating in silico at least onemathematical model representing one or more physiological systems. Themodel can include one or more adjustable parameters for each of the oneor more physiological systems. To generate a virtual patientcorresponding to an actual patient, the values of the adjustableparameters can be determined based on physiological data, including aresponse of that actual patient to administration of a therapeuticregimen.

In one aspect, the invention relates to a method for assessingmodifications to a therapy. The method comprises generating in silico aplurality of virtual patients based on data collected from a populationof previously treated patients; utilizing the virtual patients to devisea simulated therapy, determining one or more physiological parameters ofthe virtual patients in response to application of the simulated therapyto the virtual patients, modifying one or more parameters of thesimulated therapy and determining one or more physiological parametersof the virtual patients in response to application of the modifiedsimulated therapy to the virtual patients, and identifying one or morepatients, if any, that respond adversely to the modified simulatedtherapy based on the determined physiological parameters.

In one aspect, the invention relates to a method of determiningcompliance of a patient suffering from a health condition with atherapeutic protocol. The method comprises generating in silico avirtual patient representing the actual patient. A simulated therapeuticprotocol can be applied to the virtual patient, and one or morephysiological parameters of the virtual patient can be determined inresponse to the simulated therapeutic protocol. Further, one or morerespective physiological parameters of the respective actual patient canbe measured and compared with the physiological parameters determined inthe virtual patient to assess the actual patient's compliance with thetherapeutic protocol.

In some embodiments, the data collected from the population ofpreviously treated patients comprises gender, age, weight, height,ethnicity, metabolic/chemistry, complete blood count, or a combinationthereof. In some embodiments, the data collected further comprisesmedications used by the previously treated patients, past medicalhistory data, past surgical history data, or a combination thereof. Inone embodiment, the past medical history data comprises informationregarding diabetes, blood pressure/hypertension, cancer, congestiveheart failure, kidney disease, or a combination thereof.

In some embodiments, the methods and computer systems described hereinfurther comprise determining whether at least one of a plurality ofdetermined physiological parameters of at least one of a plurality ofvirtual patients is indicative of an adverse effect from a simulatedtherapy.

In some embodiments, one or more physiological parameters in a pluralityof virtual patients are determined in response to an applied therapy ata plurality of times over a predetermined time interval.

In some embodiments, the one or more physiological parameters comprisesone or more metabolic parameters. The metabolic parameters can comprise,for example, blood parameters and/or urine parameters. By way ofexample, the blood parameters can comprise any of hemoglobinconcentration and hematocrit level.

In some embodiments, each of the virtual patients comprises at least onein silico model representing a physiological system. By way of example,the model can be an erythropoiesis model, which can represent, at leastin part, the physiological system responsible for production of redblood cells.

In some embodiments, the therapy comprises any of at least onepharmacological therapy, at least one non-pharmacological therapy, or acombination thereof. In one embodiment, the at least one pharmacologicaltherapy comprises an anemia therapy. The anemia therapy can compriseadministering at least one of an erythropoiesis stimulating agent (ESA),iron, or a combination thereof. In some embodiments, the erythropoiesisstimulating agent comprises exogenous erythropoietin. In one embodiment,the one or more non-pharmacological therapies comprises a fluid therapy,a dietary therapy, an exercise therapy, or a combination thereof.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing will be apparent from the following more particulardescription of exemplary embodiments of the invention, as illustrated inthe accompanying drawings in which like reference characters refer tothe same parts throughout the different views. The drawings are notnecessarily to scale, emphasis instead being placed upon illustratingembodiments of the present invention.

FIG. 1 illustrates patient selection for 3 clinical trials from a targetpopulation.

FIG. 2A is a flow diagram of an in silico trial used to supplementactual clinical trials.

FIG. 2B is a flow diagram of one embodiment of a method according to anembodiment to perform a virtual trial.

FIG. 2C is a flow diagram of one embodiment of a computer systemaccording to the present teachings.

FIG. 2D is a flow diagram modeling a process to determine operationalgoals of a clinical trial.

FIG. 3A illustrates generation of virtual patients from a targetpopulation for one or more virtual trials.

FIG. 3B illustrates a flow diagram of generation of virtual patientsaccording to the present teachings.

FIG. 4A illustrates data collected from a population of patients shownas a user interface on a computer system.

FIG. 4B is a table of data collected from patients in a targetpopulation.

FIGS. 5A-5B illustrate hemoglobin concentration and EPO doseadministered as a function of time in a patient. In FIG. 5A, the dashedline represents the actual hemoglobin concentration of an actual patientand the solid line represents the simulated hemoglobin concentration forthe patient.

FIGS. 6A-6B illustrate hemoglobin concentration and EPO doseadministered as a function of time in a virtual patient.

FIGS. 7A-7C illustrate hemoglobin concentration and two differentsimulated EPO dose administered to a virtual patient as a function oftime.

FIGS. 8A-8C illustrate hemoglobin concentration and two differentsimulated EPO dose administered to a virtual patient as a function oftime.

FIG. 9A is a schematic illustration of the organizational diagram of ananemia (erythropoiesis) model suitable for generating virtual patientsin some embodiments of the present teachings.

FIG. 9B is a schematic illustration of the organizational diagram of ananemia (erythropoiesis) model including an iron submodel suitable forgenerating virtual patients in some embodiments of the presentteachings.

FIG. 10 is a list of equations used to represent various aspects of theaforementioned anemia (erythropoiesis) model.

FIG. 11 illustrates the interrelationship between differentphysiological systems that can be represented in a mathematical model.

FIG. 12 is a table of data from virtual patients in a virtual dialysisclinic (VDC).

FIG. 13 illustrates the erythropoiesis in humans.

FIG. 14 is a schematic view of a computer network in which someembodiments of the present invention can be implemented.

FIG. 15 is a block diagram of a computer of the network of FIG. 14 .

FIG. 16 is a schematic of a hemodialysis machine.

FIG. 17 is a perspective view of a drug delivery device of thehemodialysis machine of FIG. 16 .

DETAILED DESCRIPTION OF THE INVENTION

In one aspect, the present teachings relate generally to methods andsystems for assessing therapies in silico using one or more virtualpatients. As discussed in more detail below, each virtual patient can bea mathematical construct representing in silico one or morephysiological systems of a respective actual patient. In this manner, aplurality of virtual patients can be constructed that mimic a pluralityof actual patients. As discussed in more detail below, such virtualpatients can be employed, for example, to optimize a therapy foradministration to actual patients, reject a therapy based on observedadverse effect on one or more virtual patients, assess how modificationsin a therapy can affect actual patients, etc. In some embodiments, acollection of virtual patients generated in accordance with the presentteachings can be employed to conduct in silico clinical trials. Such insilico clinical trials can provide a number of advantages, such asinclusion of individuals (such as, pregnant women, infants) that aretypically difficult to include in clinical trials, the possibility ofmodifying a therapy and obtaining information regarding modified therapyon a much shorter time scale than is possible in conventional clinicaltrials.

Certain exemplary embodiments will now be described to provide anoverall understanding of the principles of the computer systems andmethods disclosed herein. One or more examples of these embodiments areillustrated in the accompanying drawings. Those skilled in the art willunderstand that the computer systems and methods specifically describedherein and illustrated in the accompanying drawings are non-limitingexemplary embodiments and that the scope of the present invention isdefined solely by the claims. The features illustrated or described inconnection with one exemplary embodiment may be combined with thefeatures of other embodiments. Such modifications and variations areintended to be included within the scope of the present invention.

So that the invention may more readily be understood, certain terms arefirst defined.

As used herein, the terms “virtual patient,” “avatar” and “virtualdummy” are used interchangeably and refer to a mathematical construct(model) that represents one or more physiological systems (orsubsystems) of an actual patient. The term “actual patient” as usedherein refers to a live or dead individual. As used herein, the terms“about” or “approximately” for any numerical values or ranges indicate asuitable dimensional tolerance that allows the composition, part, orcollection of elements to function for its intended purpose as describedherein. These terms indicate at most a ±5% variation about a centralvalue.

As used herein, the terms “a” or “an” refers to one or more of this kindof entity; for example, “an adverse effect” should be understood as oneor more adverse effects e.g., from a therapy (e.g., drug). Therefore,the terms “a,” “an,” “one or more” as well as “at least one” can be usedinterchangeably herein.

Virtual Dialysis Clinic (Virtual Clinical Trials)

Described herein is a Virtual Dialysis Clinic (“VDC”) or a method andsystem for performing a virtual clinical trial, in silico, for a healthcondition. For example, the methods and systems described herein canprovide a complementary approach to the current standards in clinicaltrials and health care practice (see FIG. 2A). In some cases, themethods and systems describe herein can replace conventional clinicaltrials. The VDC approach can bypass and/or address ethical, legal, andscientific challenges typically encountered in current approaches toclinical trials for determining novel therapies and diagnostics.

Clinical research regarding new therapeutics and/or diagnostics follow agenerally well-defined pathway in the form of clinical trials (Phases 1,2, 3 and 4). These trials take place in volunteers (phases 1-3), who areeither healthy (phase 1) or come from the intended target population(phases 2 and 3). Each volunteer undergoes specific pre-test screeningto assess their eligibility, using inclusion and exclusion criteria setby the study design. The most rigorous standard, a phase 3 trial, is arandomized, controlled trial (RCT). Randomization is applied to accountfor observed and unobserved baseline covariates. Further, there areclinical and statistical (e.g., power, confidence, etc.) considerationsthat inform sample size decisions (e.g., how many individuals need to beenrolled in a study to observe a positive or negative effect, how manytotal patients having a certain disease or condition should be enrolledin a clinical trial). The patient population in a phase 3 trial is, inmost cases, a highly selective group. However, some of the majorconcerns with any clinical trial are the generalizability of the studyfindings to the intended target population, population size, safety(e.g. adverse effects), and costs.

FIG. 1 schematically illustrates three conventional clinical trials. Thetarget population comprises all individuals (represented by the dots)that would be selected to enter one of the three clinical trials:clinical trial 1, clinical trial 2 or clinical trial 3. As shownschematically in the figure, each individual of the target populationcan participate in only one of the three clinical trials, thus limitingthe number of individuals available for each trial. Further, it may notbe possible to include certain individuals in the clinical trials, e.g.,due to safety and/or ethical concerns. Moreover, each clinical trial canrequire many months, or even years, to complete, and can cost hundredsand even millions of dollars. As discussed further below, the teachingsof the invention for performing clinical trials in silico address suchand other shortcomings of conventional clinical trials.

FIG. 2A represents a flow diagram of a clinical trial that issupplemented with a virtual trial in accordance with one embodiment ofthe present teachings. In Phase 1, a therapeutic (e.g., drug) is testedfor safety, dosage, and side effects in a small number of healthyvolunteers. Phase 2 examines the drug's safety and effectiveness in asmall number of volunteers (those having the disease or is in need oftreatment).

Rather than moving after phase 2 to a conventional phase 3 trial, an insilk( ) (virtual) trial (e.g., VDC) according to the present teachingscan be performed, e.g. to inform the conduct of a subsequent phase ofthe clinical trial. Such a virtual trial can include, for example,extensive high-throughput testing for safety and effectiveness employingcomputer simulations in appropriate physiological models of existingpatients. Based on information and data gathered from the virtualclinical trial, a phase 3 trial can be performed, optimized based on theresults of the virtual trial. A phase 3 clinical trial looks at theeffectiveness, side effects, and comparison with other treatments in alarge number of volunteer patients.

A virtual trial, in silico, is run in an initial step similarly to aphase 1, 2 or 3 clinical trial (see, e.g., FIG. 3A) in that a questionis typically defined to be addressed in the trial (e.g., How does betterfluid management affect erythropoiesis stimulating agents (ESA) and ironutilization?; What is the effect of iron delivery in the dialysate onESA utilization?; How does an anti-inflammatory intervention affectanemia control?). Next, mathematical models are designed and used todescribe the physiology pertinent to the trial question at hand. Arepresentative and/or random actual patient sample is drawn from alarger target population. As discussed in more detail below,mathematical techniques are employed to create patient specific virtualpatients based on the actual patients in the drawn patient population.One or more test treatment options are applied, in silico, to thevirtual patients. The responses of the virtual patients to the appliedtest treatment options can inform which option(s) can provide the mostpromising treatments for application to actual patient population. Insome embodiments, the most promising treatment options can be selectedfor further clinical (actual) testing.

For example, as shown in the flow chart depicted in FIG. 2B, in oneembodiment, the invention relates to a method for performing a clinicaltrial, in silico, (herein also referred to as a “virtual clinical trial”or for brevity a “virtual trial”) for a health condition, where themethod comprises generating in silico one or more virtual patients (step1). The virtual patients are generated based on data collected from apopulation of previously treated actual patients suffering from thehealth condition. The method further comprises applying a simulatedtherapy to each of the one or more virtual patients (step 2) anddetermining one or more physiological parameters in the virtual patientsin response to the simulated therapy (step 3).

The above method of performing a clinical trial can further compriseperforming iteratively the following steps so as to derive an optimaltherapy for the health condition. (i) modifying one or more parametersof the simulated therapy, thereby creating a modified simulated therapy(step 4); (ii) administering the modified simulated therapy to each ofthe virtual patients (step 5); and (iii) determining one or morephysiological parameters in the virtual patients in response to themodified simulated therapy (step 6).

In another embodiment, the invention relates to a computer system forperforming a clinical trial (e.g., determining efficacy of a therapy).By way of example, FIG. 2C schematically depicts a computer system 10according to an embodiment of the present teachings for determining insilico the efficacy of a therapy. The computer system 10 comprises ageneration module 11, a simulation module 21, a processor 31, a randomaccess memory (RAM) 41, an input/output interface 42, and one or moredatabases, such as database 51 and database 52. Although two databasesare depicted in this exemplary embodiment, in other embodiments a singledatabase, or more than two databases, can be employed. A bus 32, whichcan comprise a set of hardware lines, allows transferring data andinstructions among various components of the computer system. Theinput/output interface 42 allows connecting the computer system to avariety of input/output devices, such as keyboard, mouse, display,printers, speakers, etc. Further, network interface 44 allows thecomputer system 10 to communicate with other network-enabled devices.The RAM 41 in turn provides volatile storage for computer softwareinstructions and data used to implement the present teachings.

The generation module 11 generates in silico a plurality of virtualpatients based on data collected from a population of previously treatedpatients, e.g., patients suffering from a health condition. The datacollected from actual patients can be inputted into the computer system10 via one or more input/out device and can be stored in database 51and/or database 52 The generation module 11 can access the stored datato generate virtual patients as mathematical constructs based on theactual patients' data, e.g., in a manner discussed in more detail below.The generated virtual patients can also be stored in database 51 and/ordatabase 52 and can be accessed by the simulation module 21.

The simulation module 21 is configured to apply a simulated therapy tothe virtual patients to determine one or more physiological parametersof the virtual patients in response to the simulated therapy. Thedetermined physiological parameters can be stored in database 51 and/ordatabase 52. For example, virtual patients providing a mathematicalmodel of erythropoiesis can be employed to determine hemoglobinconcentration and hematocrit level in response to simulatedadministration of an EPO regimen.

The simulation module 21 can further be configured to adjust thesimulated therapy based on the one or more physiological parametersthereby creating a modified simulated therapy. The simulation module 21can apply the modified simulated therapy to the virtual patients, anddetermine one or more physiological parameters of the virtual patientsin response to the modified simulated therapy. The simulation module 21can iteratively repeat the process of modifying the therapy anddetermining one or more physiological parameters until an optimaltherapy is obtained. For example, if an initial simulated EPOadministration were to result in low hemoglobin concentration, the levelof EPO and/or its frequency of administration may be increased and thelevel of hemoglobin concentration in response to this modified therapycan be determined, to decide whether additional modifications to the EPOregimen is needed.

In some embodiments, the simulation modulation can be configured torecognize one or more physiological parameters determined in response toan applied simulated therapy that are indicative of an adverse effect.Further details regarding a computer system according to variousembodiments of the present teachings can be found further below inconnection with the description of FIG. 15 .

Referring to FIG. 2D, a modeling cycling process can be based on anumber of factors. First, goals of the modeling process can be defined(201). After empirical observations, and understanding mechanisms andrelationships of various factors (202), a mathematical model (203) canbe created. Simulations (204) using that mathematical model are thentested and validated with actual experiments (205). These results areexamined (207). If the results are not satisfactory, the model can berefined and/or redefined (206). The goals of the model can also bechanged according to the findings. However, the results can align withthe goals of the model and would then be ready for implementation (208)in a virtual trial.

Referring to FIGS. 3A-3B, a (one or more) virtual or in silico trial canbe performed on a randomized population of virtual patients. Forexample, intervention 1 can be performed on a random sample of virtualpatients, generating result 1. Intervention 2 can also be performed onthe random sample of virtual patients, generating result 2. Intervention3 can be performed on the random sample of virtual patients, generatingresult 3. Because the virtual trials are performed in silico,interventions 1, 2 and 3 can be done simultaneously or sequentially. Forexample, in some instances it is preferred to have interventions 1, 2,and 3 performed simultaneously so that results 1, 2, and 3 are known ator around the same time. In other instances, it is preferred to haveintervention 1 performed first. Intervention 2 can be then created,modified, or otherwise changed based on result 1 from intervention 1.Then, intervention 3 can also be created, modified, or otherwise changedbased on result 2 from intervention 2. For illustration purposes only,FIG. 3A illustrates 3 different interventions providing 3 results.However, it is possible to compute tens, hundreds, and even thousands ofinterventions (e.g., therapies) on the virtual patients.

FIG. 3B also depicts a diagram of performing a virtual trial accordingto the present teachings. Virtual patients (avatars) 304 can be createdand generated from data obtained from actual patients (e.g., from adialysis center 301). These virtual patients can be used in any numberof virtual clinical trials 305, where one or more interventions aretested. These interventions can be modified or changed according to theoutcomes of the virtual clinical trials. Any one of the interventionscan then be performed in actual clinics on patients based on the resultsof the virtual trials.

For example, the mean data from a population of virtual patients in avirtual dialysis clinic is shown in the table in FIG. 12 . This datacorresponds to the mean values of various parameters, such as age,weight, albumin, EPO dose, ferritin, and others.

Compared to traditional clinical trials, virtual trials according to thepresent teachings have many advantages. Virtual trials are safe becauseno actual patients are used. Also, virtual trials can be performedfaster since the in silico models can be employed to providehigh-throughput evaluations of therapeutic and/or diagnosticalternatives. In some cases, the results of virtual trials can be morereliable than those of traditional clinical trials. For example, in somevirtual trials, the absence of inclusion and/or exclusion criteria canincrease generalizability and reliability. Virtual trials can also bemore inclusive than traditional clinical trials, since classes ofpatients or individuals that may be excluded from a traditional clinicaltrial, such as pregnant women, children, and the elderly, can beincluded in a virtual trial. Further, virtual trials can be moreeconomical than conventional clinical trials, for example, largemulticenter clinical trials. In particular, no expensive trialinfrastructure is required when performing a virtual trial.

In some embodiments, a virtual trial can determine efficacy of atherapy, medical device and/or a diagnostic intervention. A virtualtrial can determine any clinical effect, any outcome, or any changesimilar to an actual clinical trial. For example, in a virtual trialaccording to the present teachings, one or more physiological parametersof the virtual patients in response to an applied simulated therapy canbe determined. The one or more physiological parameters can allow aclinician or investigator to determine efficacy of the therapy. Further,the physiological parameters may indicate little or no effect, or apositive effect, from the (simulated) therapy. The physiologicalparameters may indicate an adverse effect from the (simulated) therapy.

The method for determining efficacy of a therapy can further compriseadjusting the simulated therapy based on the determined one or morephysiological parameters of said virtual patients, applying the adjustedtherapy to the virtual patients, and determining one or morephysiological parameters in response to the adjusted simulated therapy(e.g., FIGS. 6A, 7A, and 8A). The method can further comprise repeatingthese steps so as to optimize the simulated therapy for application toone or more actual patients.

The one or more physiological parameters can be any parameter measured,determined, or otherwise known in an individual. For example, the one ormore physiological parameters can comprise one or more metabolicparameters. Metabolic parameters, for example, can comprise blood orurine parameters. By way of example, blood parameters can comprisehemoglobin, hematocrit, complete blood count (white blood cell, whiteblood cell differential, neutrophil count, lymphocyte count, monocytecount, eosinophil count), red blood cell, MCV (mean corpuscular volume),MCH (mean corpuscular hemoglobin), MCHC (mean corpuscular hemoglobinconcentration), platelet count, or combinations thereof. Other bloodparameters can comprise a comprehensive metabolic panel. A comprehensivemetabolic panel comprises serum glucose, calcium, BUN (blood ureanitrogen), creatinine, sodium, potassium, chloride, carbon dioxide,total protein, albumin, bilirubin, ALP (alkaline phosphatase), AST(aspartate aminotransferase), ALT (alanine aminotransferase), orcombinations thereof.

In a virtual trial described herein, one or more physiologicalparameters can be determined at a single point or multiple points intime. Multiple points of time can be defined over a predetermined timeinterval (e.g., 1 day, 1 month, 1 year) or can be open-ended (e.g.,indefinitely). Also, a physiological parameter can be determined at anyfrequency for a given or predetermined time interval. For example, a redblood cell count can be determined every day for as long as a certaintherapy is administered.

In a virtual clinical trial, a therapy can comprise a pharmacologicaltherapy, a non-pharmacological therapy, or a combination thereof. Forexample, a pharmacological therapy can comprise an anemia therapy, adiabetes therapy, a hypertension therapy, a cholesterol therapy,neuropharmacologic drugs, drugs modulating cardiovascular function,drugs modulating inflammation, immune function, production of bloodcells; hormones and antagonists, drugs affecting gastrointestinalfunction, chemotherapeutics of microbial diseases, and/orchemotherapeutics of neoplastic disease. Other pharmacological therapiescan include any other drug or biologic found in any drug class. Forexample, other drug classes can comprise allergy/cold/ENT therapies,analgesics, anesthetics, anti-inflammatories, antimicrobials,asthma/pulmonary therapies, cardiovascular therapies, dermatologytherapies, endocrine/metabolic therapies, gastrointestinal therapies,cancer therapies, immunology therapies, neurologic therapies, ophthalmictherapies, psychiatric therapies or rheumatologic therapies. Anon-pharmacological therapy can comprise a fluid therapy, a dietarytherapy, an exercise therapy, a medical device, diagnostic,extracorporeal therapies, radiotherapy, therapies using sound, includingultrasound, and/or electrotherapies.

In some embodiments, an anemia therapy comprises administering at leastone erythropoiesis stimulating agent (ESA), iron, drugs stimulatingendogenous erythropoietin release and synthesis (e.g. hypoxia induciblefactor (HIF) stabilizers), and/or biosimilars. The ESA can compriseexogenous erythropoietin (EPO), or synthetic EPO. In other embodiments,the EPO comprises recombinant EPO.

By way of example, FIGS. 6A through 8C illustrate graphs of hemoglobinconcentration as a function of time in virtual patients in response tosimulated administration of EPO according to different regimens. Thevirtual patients were generated based on an EPO model described inPCT/US2012/054264, entitled “System and Method of ModelingErythropoiesis and Its Management,” which is herein incorporated byreference in its entirety. Each hemoglobin concentration curve (e.g.,FIGS. 6A, 7A, and 8A) correspond to an administration of EPO, in silico,over the same time period (360 days). Each of the predicted hemoglobinconcentrations is due to a hypothetical administration of EPO. Forexample, FIGS. 6A and 6B illustrate the effect of simulatedadministration of EPO on hemoglobin concentration in a virtual patientover a number of days. The EPO dose generally tracks with hemoglobinconcentration, i.e., the higher the EPO dose, the greater the hemoglobinconcentration and the lower the EPO dose, the lower the hemoglobinconcentration.

Different simulated therapies can be administered and compared, such asthose illustrated in FIGS. 7A-7C. In FIG. 7B, a bolus-likeadministration of EPO is applied to a virtual patient for 60 days, witha break in EPO administration (i.e., no EPO administered) for about 40days. The virtual patients were generated based on an EPO modeldescribed in PCT/US2012/054264. This pattern of EPO administration isthen repeated. The corresponding hemoglobin concentration determined inthe virtual patient is shown in FIG. 7A (graph A). Graph A shows thatthe hemoglobin concentration exhibits an oscillatory behavior with thepeak of the response shifted (delayed) relative to the peak of EPOadministration and the minimum of the response shifted relative to theperiods of no EPO administration, indicating a delay in responserelative to a change in EPO administration. Graph A further shows thatthe hemoglobin concentration deviates from the desired range of 10 to 11g/dl during certain temporal periods.

In FIG. 7C, the EPO administration is generally at a lower, more evendose level than that shown in FIG. 7B. The response of the virtualpatient to the administration of EPO according to FIG. 7C is shown ingraph B in FIG. 7A. Graph B shows a more even response with lowermaximum-to-minimum variation.

Similarly, FIGS. 8A-8C also illustrate different simulated EPOadministrations (FIGS. 8B and 8C) and how they would affect a virtualpatient's hemoglobin concentration (FIG. 8A). The virtual patients weregenerated based on an EPO model described in PCT/US2012/054264.Specifically, graph A in FIG. 8A shows the temporal variation of thehemoglobin concentration in a virtual patient subjected to the EPOregimen shown in FIG. 8B and graph B in FIG. 8A shows the temporalvariation of the hemoglobin concentration in a virtual patient subjectedto the EPO regimen shown in FIG. 8C.

The data obtained by applying different EPO regimens to virtual patientscan inform selecting an optimal regimen for actual patients. Forexample, the data provided in FIGS. 7A-7C and 8A-8C discussed aboveshows that the EPO regimen depicted in FIG. 8C causes a response inhemoglobin concentration that remains within the desired range (i.e.,between 10 and 11 g/dl) over essentially the entire period of 360 days,and hence may be more suitable than the other depicted regimens foradministration to actual patients.

In some embodiments, a simulated therapy can be adjusted based on thedetermined one or more physiological parameters. Subsequently, one ormore physiological parameters can be determined from the adjustedsimulated therapy. The process can be repeated until, e.g., an optimalsimulated therapy is obtained, which results in the physiologicalparameters being in a desired range. It would be generally apparent toone of skill in the art that the physiological parameters determined inthe adjusted simulated therapy can be the same or different than thephysiological parameters determined in the simulated therapy. In certainembodiments, the one or more physiological parameters are the same inthe simulated therapy and the adjusted therapy. In other embodiments theone or more physiological parameters in the simulated therapy can bepartially or completely different than the one or more physiologicalparameters in the adjusted simulated therapy. It will also beappreciated that by virtue of adjusting a therapy, other, differentparameters (e.g., physiological parameters) can become important in thevirtual trial.

As noted above, a simulated therapy can be adjusted multiple times basedon the determined one or more physiological parameters. After eachadjusted simulated therapy, one or more physiological parameters can bedetermined. The simulated therapy can be then readjusted or refined tomaximize the outcome (e.g., efficacy) of a therapy. For example, in oneembodiment, an optimal therapy for a health condition can be ascertainedby iteratively applying a simulated therapy to a plurality of virtualpatients, determining values of one or more physiological parameters inthe virtual patients in response to the therapy, applying an adjustedsimulated therapy to the virtual patients, determining one or morephysiological parameters in response to the adjusted therapy until thephysiological parameters are within a desired range, indicating that anoptimal therapy is achieved.

In one embodiment, the invention relates to a method and system forassessing modifications to a therapy. The method comprises generating insilico a plurality of virtual patients based on data collected from apopulation of previously treated patients, utilizing said virtualpatients to devise a simulated therapy, determining one or physiologicalparameters of said virtual patients in response to application of saidsimulated therapy to said virtual patients, modifying one or moreparameters of said simulated therapy and determining one or morephysiological parameters of said virtual patients in response toapplication of said modified simulated therapy to said virtual patients,and identifying one or more patients, if any, that respond adversely tosaid modified simulated therapy based on said determined physiologicalparameters.

In the above method, the identification of a therapy that leads toadverse effects in one or more patients can be useful in eliminatingthat therapy from a set of therapies that can be recommended foradministration to actual patients.

In one embodiment, the invention relates to a method of determining theeffects of non-compliance with a therapeutic protocol in a patientsuffering from a health condition. The method comprises generating insilico a virtual patient representing said patient, wherein said virtualpatient models said health condition. The method further comprisesapplying a simulated therapeutic protocol that takes non-compliance intoaccount to said virtual patient, determining one or more physiologicalparameters of said virtual patient in response to said simulatedtherapeutic protocol, measuring corresponding one or more physiologicalparameters of a patient, and comparing said measured one or morephysiological parameters with said determined one or more physiologicalparameters to assess the effects of non-compliance with said therapeuticprotocol.

The results of the simulations (e.g., simulated therapy) can be used tomake informed decisions about the actual roll-out of specific treatmentoptions. Multiple computer simulations of treatment options in apopulation can generate extensive amounts of data that would otherwisetake months to years to collect in an actual clinical trial. Also,results can be extrapolated to the target population. The most promisingtreatment options can be selected for further clinical testing oradministration to actual patient without further testing.

Virtual Patients

As discussed above, the methods according to the present teachings relyon the generation of patient-specific virtual patients. As discussed inmore detail below, each virtual patient is generated so as to mimic atleast one physiological system (or subsystem) of an actual patient. Inother words, in many embodiments, there is a one-to-one correspondencebetween a virtual patient and an actual patient based on which thevirtual patient is generated. As discussed in more detail below, avirtual patient can be represented by a mathematical model whoseadjustable parameters are determined based on data collected from anactual patient. By way of example, the data can be measured hemoglobinconcentration in an actual patient suffering from diabetes in responseto administration of a regimen of EPO to that patient.

As noted above, a virtual patient is generated, in silico, based on datacollected from an actual patient (see, e.g., Table 1 below). A pluralityof virtual patients can be created based on data collected from apopulation of actual patients. For example, FIGS. 4A-4B represent anaugmented data set of the VDC, including previously collected data fromactual patients (e.g., a dummy patient ID, gender, age, race, BMI,diabetes status) and parameter values determined using advancedmathematical techniques (RBC lifespan, BM (bone marrow) reaction,endogenous EPO concentration and EPO half-life). Each actual patient ina target population can be represented as a mathematical construct(e.g., as a virtual patient). The collected data can represent at leastone measured biological response of the previously treated patients to apreviously administered therapeutic regimen.

TABLE 1 Biological characteristics. Biological key characteristics“Normal” TBV 3567 mL Stem cells 5.92*10{circumflex over ( )}6 RBClifespan 92 days EPO half-time 4.8 h End. Edo 5.3 mU/mL Bone marrowreaction Normal

Each virtual patient can mathematically represent one or morephysiological systems of an actual patient. Previously collected datafrom the actual patient allows for the construct of such a virtualpatient. For example, the collected data can represent at least onemeasured biological response in the actual patient to a previouslyadministered therapeutic regimen. The data collected from a populationof previously treated actual patients can further comprise gender, age,weight, height, ethnicity, metabolic parameters, chemistry parameters,complete blood count, or combinations thereof. Any clinically relevantdata or patient information can be used in the construct of a virtualpatient. For example, medications taken by the actual patient, pastmedical history, past surgical history, family history, or combinationsthereof can also comprise the data collected from the population ofactual patients. Past medical information can include informationregarding diabetes, blood pressure/hypertension, cancer, congestiveheart failure, kidney disease. Past medical information can also includeany information found in the International Statistical Classification ofDiseases and Related Health Problems (ICD), such as the ICD-9,including, diseases of the blood and blood-forming organs and certaindisorders involving the immune mechanism; endocrine, nutritional andmetabolic diseases; mental and behavioral disorders; diseases of thenervous system; diseases of the eye and adnexa; diseases of the ear andmastoid process; diseases of the circulatory system; diseases of therespiratory system; diseases of the digestive system; diseases of theskin and subcutaneous tissue; diseases of the musculoskeletal system andconnective tissue; diseases of the genitourinary system; pregnancy,childbirth and the puerperium; certain conditions originating in theperinatal period; congenital malformations, deformations and chromosomalabnormalities; injury, poisoning and certain other consequences ofexternal causes; external causes of morbidity; infectious and parasiticdiseases. Collected data from a population of patients can furthercomprise any information found in their medical records or patientcharts.

As previously described, in one embodiment, one or more mathematicalmodels can be used to create a virtual patient. Mathematical models canbe used to represent any biological or physiological system, componentor subcomponent. For example, any one of the following biologicalsystems can be represented mathematically: circulatory system (lungswith heart, blood and blood vessels), integumentary system (skin, hair,fat, and nails), skeletal system (structural support and protection withbones, cartilage, ligaments and tendons), reproductive system (the sexorgans, such as ovaries, fallopian tubes, uterus, vagina, mammaryglands, testes, vas deferens, seminal vesicles and prostate), digestivesystem (digestion and processing food with salivary glands, esophagus,stomach, liver, gallbladder, pancreas, intestines, rectum and anus),urinary system (kidneys, ureters, bladder and urethra involved in fluidbalance, electrolyte balance and excretion of urine), respiratory system(the organs used for breathing, the pharynx, larynx, bronchi, lungs anddiaphragm), endocrine system (hormones made by endocrine glands such asthe hypothalamus, pituitary gland, pineal body or pineal gland, thyroid,parathyroid and adrenals, i.e., adrenal glands), lymphatic system(structures involved in the transfer of lymph between tissues and theblood stream; includes the lymph and the nodes and vessels, includesfunctions including immune responses and development of antibodies),muscular system (manipulation of the environment, provides locomotion,maintains posture, and produces heat; includes skeletal muscles, smoothmuscles and cardiac muscle), and nervous system (collecting,transferring and processing information with brain, spinal cord andperipheral nervous system).

Specific models (modules) can be created for more specific physiologicalsystems, such as, for example, bone mineral metabolism, acid-basemetabolism, cardiopulmonary system and tissue oxygenation, and sodiumand other mineral homeostasis.

Other physiological systems to which a model (or module) can be createdinclude cardiovascular system; endocrine system; models related to bloodcell physiology and immunology; inflammation; hepatic andgastrointestinal function; renal function; and/or musculo-skeletalfunction.

For example, a virtual patient can be constructed using an anemia modeldescribed herein. The model can be fitted to actual patient data (e.g.,hemoglobin concentration) to obtain values for certain parameters of themodel. In some embodiments, other parameters, such as blood volume andthe number of stem cells committed to the erythroid lineage can beestimated based on gender, weight, and height of the patient. Otherparameters can be adapted for each individual patient, such as, redblood cell life span, endogenous EPO production, half-life ofepoetin-alfa, bone marrow reaction to EPO, rate of apoptosis of CFU-Ecells, and maturation velocity of reticulocytes.

By way of example, measured data indicating the hemoglobin concentrationof a diabetic actual patient in response to administration of an EPOregimen to that patient can be used to determine values for certainparameters of a virtual patient representing that actual patient. Forexample, such data can be utilized to determine one or more adjustableparameters of a mathematical model for Erythropoiesis. For example, theparameters can be adjusted so as to obtain a simulated response (e.g., aresponse characterized by hemoglobin concentration) to simulatedadministration of the shown EPO regimen that is sufficiently close tothe actual response. In some cases, a merit function can be employed toquantify the difference between the simulated response and the actualresponse. The parameters can be adjusted until the merit functionindicates an acceptable difference between the simulated and the actualresponses.

By way of further illustration, FIGS. 5A and 5B, represent,respectively, measured hemoglobin concentration over a period of timeand EPO administered to that patient over the same period of time. Thepatient was a female individual, 163 cm in height, and about 66.9-68.2kg in weight during the course of treatment. From this patientinformation, it was determined that the average red blood cell lifespanwas approximately 63.2 days, endogenous EPO production was greater than20 mU/ml (milli-units per milliliter), the bone marrow was severelydepressed and the half-life of Epoetin alfa was about 4.4 hours. Thisdata was utilized for determine values for certain parameters of amathematical model for generating red blood cells (by adjusting theparameters to obtain a desired fit to the data).

As noted above, for each patient, a patient specific virtual patient oravatar is created. A database of virtual patients can be created from adatabase of actual patients. From the database of virtual patients, oneor more virtual patients can be selected for each virtual clinicaltrial. In some cases, specific virtual patients can be selected todetermine the answers to questions that would be difficult or impossibleto answer in an actual clinical trial. As noted above, individuals thatwould not be typically included in a clinical trial can be easilyincluded in a virtual trial. Virtual patients who are, for example, veryold, very young, and/or pregnant can be included in a virtual trial.

For example, an anemia therapy (erythropoiesis) model can be used tocreate a virtual patient to study medicines and therapies effecting redblood cell production, hemoglobin, hematocrit, or combinations thereof.Each of the essential aspects of erythropoiesis can be represented ineach virtual patient, described herein. The virtual patients can then besubjected to one or more anemia treatment schemes. Other mathematicalmodels can be used to represent other physiological systems. Forexample, as shown in FIG. 10 , mathematical models can be used torepresent anemia (erythropoiesis) fluid status, inflammation, and/oriron metabolism. Additional details regarding the erythropoiesis modelis described below.

Further details regarding an erythropoiesis model suitable forgenerating a virtual patient according to the present teachings formodeling red blood cell production in actual patients can be found inU.S. Published Patent Application No. 2014/0128791, entitled “System andMethod of Modeling Erythropoiesis Including Iron Homeostasis,” which isherein incorporated by reference in its entirety.

Referring to the scheme illustrated in FIG. 11 , in some embodiments,various modules representing physiological systems of an actual patientcan interact in silico to obtain a better understanding of the responsefor a patient to various therapeutic protocols. For example, thedepicted model can allow answering questions such as. How does betterfluid management affect ESA and iron utilization? What is the effect ofiron delivery in the dialysate on ESA utilitization? How does ananti-inflammatory intervention affect anemia control?

Additional modules can comprise any module for a physiological system orpathway. For example, there can be a module for bone mineral metabolism,acid-base metabolism, cardio pulmonary system, renal and/or sodium modelin, e.g., chronic kidney disease.

Anemia Therapy (Erythropoiesis) Model Red Blood Cell Physiology

Red blood cells (erythrocytes) are essential for the distribution ofoxygen through the body to organs and tissues. They take up oxygen inthe lungs and deliver it to tissues while squeezing through thecapillaries. To fulfill this task properly, they are highly specialized.For instance, being shaped like biconcave disks optimizes the oxygenexchange. Furthermore, the nuclei, organelles, and mitochondria havebeen expelled in earlier developmental stages, and therefore more spacefor hemoglobin, the molecule which oxygen binds to, is available.Erythrocytes are very deformable and can therefore pass capillaries halftheir diameter. During microcirculation, they have to withstand highshear stresses, rapid elongation, folding, and deformation. Over time,the cell membrane is damaged by these extraordinary stresses. Because ofthe lack of nuclei and organelles, red blood cells cannot divide orrepair their cell membranes. Senescent erythrocytes lose theirflexibility due to their fragmented membranes. These stiff cells coulddo harm to small capillaries or even clog them. To avoid this potentialharm, old erythrocytes are recognized by phagocytes and destroyed. Thisphagocytosis mainly takes place in the spleen and cells of thereticulo-endothelial system (RES). See Jandl, J. H., Blood: Textbook ofHematology, 2nd Ed. Little, Brown and Company, 1996 (hereinafter“Jandl”).

To compensate for phagocytosis of senescent red blood cells, it isnecessary to build new erythrocytes continuously. The maturation ofundifferentiated stem cells to mature erythrocytes is callederythropoiesis and takes place in the bone marrow. Erythropoiesis notonly has to compensate for the continuous loss of old erythrocytes, butalso for the additional loss of cells due to random breakdown, as wellas due to internal and external bleeding. Furthermore, the number of redblood cells has to be adjusted to varying environmental conditions, asfor instance a transition from low to high altitudes or vice versa, byincreasing the rate of erythropoiesis, or, conversely, by neocytolysis,a process believed to be wherein macrophages start to phagocytose youngerythrocytes (neocytes).

During the process of erythropoiesis, the cell population undergoes aseries of proliferations and differentiations (see FIG. 13 ). Startingfrom multipotential stem cells, erythroid cells mature to BFU-Es(earliest stage of erythroid committed cells), CFU-Es, different stagesof erythroblasts, and finally reticulocytes. The reticulocytes arereleased from the bone marrow into blood and mature within 1-2 days toerythrocytes (see FIG. 9A).

The primary control of erythropoiesis is governed by the hormoneerythropoietin (EPO). EPO is released into the blood stream by thekidneys based on a negative feedback mechanism that reacts to thepartial pressure of oxygen in blood. The concentration of EPO affectsthe number of circulating red blood cells by determining the number ofcells that mature into erythrocytes, either by recruitment or bypreventing apoptosis (programmed cell death), and by affecting thevelocity of maturing of progenitor and precursor cells. Thus,disturbances in oxygen delivery can be adjusted for by an adaptiveresetting of the rate of erythropoiesis. Additionally, as alreadymentioned above, there exists a physiological process which affects theselective degradation of young erythrocytes in situations of red cellexcess, called neocytolysis. Neocytolysis seems to be triggered by adrop in the EPO level. See Rice, L., and Alfrey, C. P., CellularPhysiology and Biochemistry, 2005, Vol. 15, pp. 245-250 (hereinafter“Rice 2005”); Rice, L., Alfrey, C. P., Driscoll, T., Whitley, C. E.,Hachey, D. L., and Suki, W., Amer. J. Kidney Diseases, 1999, Vol. 33,pp. 59-62 (hereinafter “Rice 1999”); and Rice, L., W. Ruiz, W.,Driscoll, T., Whitley, C. E., Tapia, R., Hachey, D. L., Conzales, G. F.,and Alfrey, C. P., Annals Internal Medicine. 2001, Vol. 134, pp. 652-656(hereinafter “Rice 2001”).

Another critical factor for effective erythropoiesis is the availabilityof iron which is indispensable for hemoglobin synthesis. If the body isnot able to provide sufficient iron for this process, then ineffectiveerythropoiesis will result. See Finch, S., Haskins, D., and Finch, C.A., The Journal of Clinical Investigation. 1950, Vol. 29, pp. 1078-1086;and Lichtman, M. A., E. Beutler, E., Kipps, T. J., Seligsohn, U.,Kaushansky, K., and Prchal, J. T. (editors), Williams Hematology, 7thedition, New York, McGraw-Hill, 2005 (hereinafter “WilliamsHematology”). In normal subjects, the total iron content of the bodystays within narrow limits (iron overload is toxic). Once an atom ofiron enters the body it is conserved with remarkable efficiency and canremain in the body for more than ten years. Iron is lost via loss ofcells (especially epithelial cells), bleeding and loss of very smallamounts via urine and sweat. The balance of iron content is achieved byabsorption and not by control of excretion. If the plasma concentrationof iron is too low, then the level of the hormone hepcidin is decreased.The consequence of a lower hepcidin level is that more iron is taken upvia the duodenum and more iron is released from macrophages and from thestores. See Crichton, R., Iron Metabolism. From Molecular Mechanisms toClinical Consequences, New York, J. Wiley, 2009 (hereinafter“Crichton”); Fleming, M. D., J. Amer. Soc. Hematology, 2008, pp. 151-158(hereinafter “Fleming”). Patients suffering from inflammation, such asdialysis patients, typically have higher hepcidin levels. Increasingiron availability in inflamed dialysis patients can be achieved by anincrease of parenteral iron by increasing dose, frequency, or both, andby reducing inflammation by diagnosis and treatment of sources ofinflammation, e.g., barrier breakdown (i.e., skin, periodontal disease,intestinal congestion), pulmonary or urinary tract infection, thrombosedfistulas or catheter, and by subsequent specific therapy, e.g.,antibiotics, catheter removal, aseptic techniques when manipulatingin-dwelling catheters, and surgical debridement of skin ulcers.

Since individual cells in the various cell populations which have to beconsidered have to be distinguished according to their age,age-structured population models are needed in order to describe thedevelopment of the cell populations. Besides these age-structuredpopulation models, the model of this invention includes a feedback loopincluding erythropoietin.

In one embodiment, a method of adjusting a patient's hematocrit and/orhemoglobin concentration to a desired range at a predetermined time withan erythropoiesis stimulating agent (ESA) regimen includes obtainingpatient parameters required for input into a model for predicting thepatient's hematocrit and/or hemoglobin concentration at a predeterminedtime with a selected ESA administration regimen, the model includingiron homeostasis, and employing the patient parameters and an initiallyselected ESA administration regimen in the model to predict thepatient's hematocrit and/or hemoglobin concentration at thepredetermined time with the initially selected ESA administrationregimen. Examples of ESAs are provided in Table 2 (adapted fromPhurrough S, Jacques L, Ciccanti M, Turner T, Koller E, Feinglass S:“Proposed Coverage Decision Memorandum for the Use of ErythropoiesisStimulating Agents in Cancer and Related Neoplastic Conditions”; Centersfor Medicare and Medicaid Services; Administrative File: CAG #000383N;May 14, 2007). See also Pfeffer, M. A., Burdmann, E. A., Chen, C-Y.,Cooper, M. E., de Zeeuw, D., Eckardt, K-U., Feyzi, J. M., Ivanovich, P.,Kewalramani, R., Levey, A. S., Lewis, E. F., McGill, J. B., McMurray, J.J. V., Parfrey, P., Parving, H-H., Remuzzi, G., Singh, A. K., Solomon,S. D., Toto, R., The New England Journal of Medicine, 2009, 361(21), pp.2019-2032; and Singh, A. K., Szczech, L., Tang, K. L., Barnhart, H.,Sapp, S., Wolfson, M., Reddan, D., The New England Journal of Medicine,2006, 355(20), pp. 2085-2098 (hereinafter “Singh et al.”). Note thatunlike other ESAs listed in Table 2, Peginesatide is not a biologicallyderived EPO; it is a synthetic peptide that stimulates EPO receptors.

TABLE 2 Erythropoiesis Stimulating Agents (EPO = erythropoietin).Compound Drug Names Manufacturer EPO α Epogen Amgen EPO α Procrit AmgenEPO α (w/o serum albumin) Eprex, Epypo J&J subsidiary (Otho Epopen,Epoxitin, Biologics) Globuren EPO β (Neo)Recormon Roche EPO β ErantinBoehringer Mannheim (Spain), Roche (Spain) EPO β Epoch Chugai EPO δ inhuman cell lines Dynepo Gene Aventis Transkaryotic Activated EPOTherapies EPO Ω Epomax, Hemax, Baxter Hemax-Eritron Modified EPO αDarbepoietin Aranesp Amgen Modified EPO α Darbepoietin Nespo AmgenModified EPO β Continuous Mircera Roche EPO Receptor Activator(Pegylation) Peginesatide Omontys Affymax

Optionally, if the patient's hematocrit and/or hemoglobin concentrationis not predicted by the model to be in the desired range at thepredetermined time with the initially selected ESA administrationregimen, the method includes employing the model with one or moredifferent ESA administration regimens until the model predicts that thepatient's hematocrit and/or hemoglobin concentration will be in thedesired range at the predetermined time. The method then includesadministering ESA to the patient with an ESA administration regimenpredicted to adjust the patient's hematocrit and/or hemoglobinconcentration to a value within the desired range at the predeterminedtime. The ESA administration regimen can include iron administration,and the model includes iron homeostasis with iron administration, andmodeling the influence of iron homeostasis on erythropoiesis. Thepatient parameters can include, for example, the starting hematocritand/or hemoglobin concentration in the patient's blood, the total bloodvolume of the patient, the lifespan of red blood cells (RBCs) of thepatient, the mean corpuscular volume of the RBCs, the rate ofneocytolysis in the patient's blood, the distribution of cell hemoglobincontent in the patient's blood, the distribution of red blood cell size,iron storage level of the patient, transferrin saturation (TSAT), serumtransferrin receptors, red cell zinc protoporphyrin level, hypochromicred blood cells, hypochromic reticulocytes, ferritn, endogenous EPOlevels, C-reactive protein (CRP) levels, interleukin-6 (IL-6) levels,and the hepcidin levels of the patient.

The starting hematocrit and/or hemoglobin concentration in the patient'sblood can be obtained from routine laboratory measurements known in theart. The total blood volume (BV) of the patient can be estimated asdescribed further below, or measured by use of radio-labeling red bloodcells with chromium-51 to estimate red blood cell volume (RCV) and usingthe formula

BV=RCV/(0.9*Hctv)

where Hctv is the venous hematocrit, obtained from routine laboratorymeasurements known in the art. See Albert S. N., Blood volumemeasurement, In: Nuclear Medicine In Vitro. 2 ed. Philadelphia: JBLippincott Co., 1983; Bernard P. J., Nouv Rev Fr Hematol 1994, 36(2),pp. 155-157; and International Committee for Standardization inHaematology: Recommended methods for measurement of red-cell and plasmavolume, J Nucl Med 1980, 21(8), pp. 793-800. The lifespan of RBCs of thepatient can be estimated from endogenous alveolar carbon monoxideconcentrations. See Strocchi A., Schwartz S., Ellefson M., Engel R. R.,Medina A., Levitt M. D., J Lab Clin Med 1992, 120(3), pp. 392-399. Themean corpuscular volume can be obtained from routine laboratorymeasurements known in the art. The rate of neocytolysis in the patient'sblood can be estimated from correlations with reduced expression of CD44(homing-associated cell adhesion molecule) and CD71 (transferrinreceptor). See Chang C. C., Chen Y., Modi K., Awar O., Alfrey C., RiceL., J Investig Med 2009, 57(5), pp. 650-654.

The models of this invention can track the patient's predictedhematocrit and/or hemoglobin concentration over time, such as betweenabout 5 days and about 200 days of the ESA administration regimen. Thepredetermined time can be any future time after an ESA administrationregimen is selected and the predicted regimen is initiated. In someembodiments, the patient undergoes a medical procedure prior, during, orafter the ESA administration regimen, such as blood donation, surgery,and dialysis, or any combination thereof. For dialysis patients, thedesired hematocrit is typically in the range of between about 28 percentand about 36 percent and the desired hemoglobin concentration istypically in a range of between about 9.5 g/dL and about 12 g/dL. SeeDriieke, T. B., Locatelli, F., Clyne, N., Eckardt, K-U., Macdougall, I.C., Tsakiris, D., Burger, H-U., Scherhad, A., The New England Journal ofMedicine, 2006, 355(20), pp. 2071-2084; Parfrey, P. S., 2006, AmericanJournal of Kidney Diseases, 47(1), pp. 171-173; Strippoli, G. F. M.,Craig, J. C., Manno, C., Schena, F. P., Journal of the American Societyof Nephrology, 2004, 15, pp. 3154-3165 (hereinafter “Strippoli et al.”);Volkova, N., Arab, L., Evidence-Based Systematic Literature Review ofHemoglobin/Hematocrit and All-Cause Mortality in Dialysis Patients,2006, 47(1), pp. 24-36. For elective orthopaedic surgery patients, thedesired hemoglobin concentration for males and females is typicallygreater than or equal to 13 g/dL, and 12 g/dL, respectively. SeeGoodnough L. T., Maniatis A., Earnshaw P., Benoni G., Beris P., BisbeE., Fergusson D. A., Gombotz H., Habler O., Monk T. G., Ozier Y,Slappendel R., and Szpalski M., British Journal of Anaesthesia 106 (1)pp. 13-22 (2011).

In another embodiment, a method of determining a patient's hematocritand/or hemoglobin concentration within a desired range at apredetermined time with an erythropoiesis stimulating agent (ESA)regimen includes obtaining patient parameters required for input into amodel for predicting the patient's hematocrit and/or hemoglobinconcentration at a predetermined time with a selected ESA administrationregimen, the model including iron homeostasis, and employing the patientparameters and an initially selected EPO administration regimen in themodel to predict the patient's hematocrit and/or hemoglobinconcentration at the predetermined time with the initially selected ESAadministration regimen. Optionally, if the patient's hematocrit and/orhemoglobin concentration is not predicted by the model to be in thedesired range at the predetermined time, or a different ESAadministration regimen is desired due to other considerations, themethod includes employing the model with one or more different ESAadministration regimens until the model predicts that the patient'shematocrit and/or hemoglobin concentration will be in the desired rangeat the predetermined time. The method then can include administering ESAto the patient with the ESA administration regimen predicted to adjustthe patient's hematocrit and/or hemoglobin concentration to a valuewithin the desired range at the predetermined time.

In yet another embodiment, a hemodialysis system includes a hemodialysismachine including a blood pump, a modular drug delivery devicecomprising at least one erythropoiesis stimulating agent (ESA) pump, atleast one iron replacement product pump, at least one ESA vial holder,and at least one iron replacement product vial holder, and a housing towhich the blood pump and the modular drug delivery device are secured.

Referring to FIGS. 1-2 (reproduced herein as FIGS. 16-17 ) of U.S. Pat.No. 8,382,696 B2 issued to Beiriger et al. on Feb. 26, 2013, which isherein incorporated by reference in its entirely, a hemodialysis system100 includes a hemodialysis machine 101 that has a drug delivery system102. The drug delivery system 102 includes a modular drug deliverydevice 103 and a disposable drug administration fluid line set 109 thatis connected to the drug delivery device 103. A drug delivery line 104of the drug administration fluid line set 109 is fluidly connected to ablood circuit of the hemodialysis system 100. The blood circuit of thehemodialysis system 100 includes, among other things, a blood line setcomprising multiple blood lines 105, a drip chamber 106, and a dialyzer110. A blood pump (e.g., a peristaltic pump) 108 is configured to pumpblood through the blood circuit during treatment. The hemodialysissystem 100 also includes a dialysate circuit and various othercomponents that, for the sake of simplicity, are not described indetail. During hemodialysis treatment, blood is drawn from the patientand, after passing through the drip chamber 106, is pumped through thedialyzer 110 where toxins are removed from the blood and collected indialysate passing through the dialyzer. The cleansed blood is thenreturned to the patient, and the dialysate including the toxins(referred to as “spent dialysate”) is disposed of or recycled andreused. As discussed in greater detail below, during the hemodialysistreatment, drugs (such as erythropoiesis stimulating agents, e.g.,Epogen®, and iron replacement products, e.g., Venofer®) are alsodelivered to the drip chamber 106 using the drug delivery system 102.The drugs mix with the patient's blood within the drip chamber 106 andare then delivered to the patient along with the patient's blood.

The modular drug delivery device 103 includes a drug vial holder 112that defines four channels 114. Each of the channels 114 is designed tohold captive a drug vial 116, 118. The channels 114 can, for example, berecesses that are formed within the drug vial holder 112 and that aresized and shaped to receive only the caps and narrow neck portions ofthe vials 116, 118 such that the larger body portions of the vials sitabove the holder 112. In the illustrated implementation, the vial 116furthest to the left contains Venofer® and the three vials 118 to theright of the Venofer® vial 116 contain Epogen®. Venofer® (iron sucroseinjection, USP) is a sterile, aqueous complex of polynuclear iron(III)-hydroxide in sucrose that is manufactured by American Regent, Inc.Venofer® is indicated in the treatment of iron deficiency anemia inpatients undergoing chronic hemodialysis who are receiving supplementalerythropoietin therapy. Epogen® is a drug that stimulates the productionof red blood cells and is also commonly used in dialysis patients.Epogen® is manufactured by Amgen, Inc.

The disposable drug administration fluid line set 109 is fluidlyconnected to each of the vials 116, 118. The drug administration fluidline set 109 includes four drug vial spikes 120 that connect to thevials 116, 118 in a manner to allow the drugs within the vials (e.g.,the Venofer® and Epogen®) to flow into feeder lines 122 via the drugvial spikes 120. Each of the feeder lines 122 is attached to aT-connector 124. The T-connectors 124 and associated tubing segments 126connect the feeder lines 124 to the drug delivery line 104. The drugvial spikes 120 can be formed of one or more relatively rigid medicalgrade plastics, such as polycarbonate or alphamethylstyrene (AMS), andthe various fluid lines can be formed of a more flexible medical gradeplastic, such as polyvinylchloride (PVC).

Each of the feeder lines 122, passes through (e.g., is threaded through)a bubble detector 128. The bubble detectors 128 are capable of detectingair bubbles within the feeder lines 122. As a result, each of the bubbledetectors 128 can determine whether its associated drug vial 116, 118 isempty during treatment, because air is drawn from the vial 116, 118 intothe feeder line 122 when the vial is empty. In some implementations, thebubble detectors 122 are optical detectors. The OPB 350 bubble detectormade by Optek can, for example, be used. Other types of opticaldetectors can alternatively or additionally be used. Similarly, othertypes of sensors, such as sensors utilizing ultrasound technology can beused as the bubble detectors. Examples of such sensors include theAD8/AD9 Integral Ultrasonic Air-InLine, Air Bubble Detector and theBD8/BD9 Integral Ultrasonic Air Bubble, Air-In-Line & Liquid LevelDetection Sensors (manufactured by Introtek International (Edgewood,N.Y.)). In some implementations, the bubble detector 128 includes asensor that, in addition to sensing the presence of an air bubble withinits associated feeder line 122, can sense the presence of the feederline itself.

Downstream of the bubble detectors 128, the feeder lines 122 passthrough (e.g., are threaded through) occluders 130. Each of theoccluders 130 can be used to crimp the portion of the feeder line 122disposed therein to prevent fluid from passing through the feeder line122. In some implementations, the occluders 130 are solenoid based rams.Alternatively or additionally, other types of automated occluders can beused. The occluders 130 can be collectively operated in a manner suchthat only one feeder line 122 is unclamped at any particular time.

The drug delivery line 104 to which each of the feeder lines 122 isfluidly connected passes through (e.g., is threaded through) aperistaltic drug pump 132. The drug pump 132 includes multiple rollersthat compress the drug delivery line 104 in a manner to create a“pillow” of fluid (i.e., a “pillow” of air or liquid) that is pinchedbetween two points of the drug delivery line 104 that are compressed bythe pump rollers. The rollers are arranged around a circumference of arotatable frame. As the frame is rotated, the rollers force the “pillow”of fluid through the drug delivery line 104 toward the drip chamber 106.When the pump 132 is being operated and one of the occluders 130 is open(i.e., not clamping its associated feeder line 122), vacuum pressure isapplied to the drug vial 116, 118 that is connected to the feeder line122 associated with the open occluder 130. In certain cases, the initialpressure in the drug vial 116, 118 is equal to the ambient pressure andwhen all of the drug has been delivered, the ending pressure within thevial is about −10 psi. In other words, the pressure within the drug vial116, 118 progresses from ambient to −10 psi as the drug is delivered.The pump 132 is configured to generate a vacuum pressure within the drugdelivery line 104 and feeder line 122 that exceeds the competing vacuumwithin the drug vial 116, 118. As a result, the drug is drawn from thevial 116, 118, through the drug vial spike 120, through the feeder line122, and into the drug delivery line 104.

In some implementations, each channel 114 of the drug vial holder 112includes a sensor to sense the presence of a vial or drug container. Incertain implementations, each drug channel 114 includes a system whichidentifies the drug vial installed. The drug vial identification systemcan, for example, include a bar code reader that reads bar codes on thevials. Different types of sensors can alternatively or additionally beused. In some implementations, for example, the vial identificationsystem uses RFTD technology. Other examples of suitable sensors includecolor sensors for sensing the color of color coded drug vials and/or forsensing the color of the drug within the vial, photo sensors (e.g.,cameras) that are equipped with text recognition software to read texton the drug vial, capacitive sensors that permit different size vials tobe detected, load cells or scales that detect the mass of the vial, andconductivity or electrical impedance sensors that can be used todetermine the type of drug within the vial.

The drug delivery device 103 also includes a control unit (e.g., amicroprocessor) that can power the various components of the drugdelivery device 103. The control unit can receive signals from and sendsignals to the various components of the drug delivery device 103,including, but not limited to, the bubble detectors 128, the occluders130, the drug pump 132, the drug vial ID sensors, and other sensorsalong the drug lines. The control unit can control the variouscomponents of the drug delivery device 103 based on information receivedfrom these components. For example, the control unit can control theoccluders 130 to ensure that only one of the occluders 130 is open at atime. This helps to ensure that drug is pulled from only one of thevials 116, 118 at a time during treatment. The control unit can alsodetermine the volume of drug delivered based on operation data of thedrug pump 132 and can control the occluders 130 based on the drug volumedetermined to have been delivered. For example, upon determining thatthe prescribed volume of the drug has been delivered, the control unitcan close the occluder 130 associated with that drug vial 116, 118 andopen the occluder 130 associated with the next drug to be delivered.

The control unit can also control the timing with which the variousoccluders 130 are opened and closed. For example, after the fullcontents of a vial have been evacuated, air will be sucked into thefeeder line 122 associated with that vial. As the air passes through thefeeder line 122, the bubble detector 128 will detect the air andtransmit a signal to the control unit indicating that the vial is empty.In response, the control unit can close the occluder 130 associated withthe empty vial and open the occluder 130 associated with the vialcontaining the next drug to be delivered. Upon receiving informationfrom the bubble detectors 128 indicating that all of the vials have beenemptied, the control unit can turn off the drug pump 132.

The control unit can also control certain components of the drugdelivery device 103 based on signals received from the drug vial IDsensors, which indicate the presence of a vial and/or the identity ofthe vial contents. Such an arrangement can help to ensure that thecorrect vials (e.g., the correct number of vials and the vialscontaining the correct contents) are used for the treatment. Uponreceiving signals from the drug vial ID sensors that do not match theinputted treatment information, for example, an alarm (e.g., an audibleand/or visual alarm) can be activated. Alternatively or additionally,the drug delivery device 103 can be configured so that treatment cannotbe initiated until the sensors detect the correct combination of vials.

The control unit can also be configured to adjust the patient'sundesired hematocrit and/or hemoglobin concentration to a value within adesired range at a predetermined time by i) employing a model forpredicting the patient's hematocrit and/or hemoglobin concentration at apredetermined time, the model including iron homeostasis, and by ii)operating the blood pump, the at least one ESA pump, and the at leastone iron replacement product pump to deliver blood, ESA, and ironreplacement product to the drip chamber via the at least one of theblood lines connected to the drip chamber and the at least two of thedrug lines connected to the drip chamber, according to an ESAadministration regimen predicted to adjust the patient's hematocritand/or hemoglobin concentration to a value within the desired range atthe predetermined time. Delivery of ESA and iron replacement product tothe drip chamber can be simultaneous, consecutive, or independent of oneanother.

The drug delivery device 103 (e.g., the control unit of the drugdelivery device 103) is configured to sense if the blood pump 108 of thedialysis machine 101 is running and to pause drug delivery if the bloodpump 108 is stopped. This technique prevents ‘pooling’ of the delivereddrug in the drip chamber 106 during treatment.

The drug delivery device 103 further includes a user interface 134 thatis connected to the control unit. The user interface 134 includesright/left arrow keys that allow the user to navigate through displaysassociated with the vials 116, 118. The user interface 134 also includesup/down arrow keys that enable the user to set the desired dosage foreach of the vials 116, 118. In addition, the user interface 134 includesstart and stop keys that allow the user to start and stop the drugdelivery device 103.

Any of various other types of user interfaces can alternatively oradditionally be used. In some implementations, the drug delivery deviceincludes a user interface that allows the user to select a drug toinfuse from a menu. In certain implementations, the user may confirmthat the drug identified by the drug vial ID sensor is correct and/ormake appropriate adjustments. The user interface can be used to inputand/or monitor various different treatment parameters. Examples of suchparameters include drug dosage, drug delivery rate, amount of drugdelivered, status of the drug delivery for each drug channel, time,percent complete, percent remaining, time remaining, time delivered,date, patient ID, patient name, alarms, alerts, etc. Such userinterfaces can include a color graphical display. In certainimplementations, for example, the user interface is color codedaccording to drug, dosing, or status of drug delivery (e.g., done,running, ready, etc.).

The drug delivery device 103 also includes an alarm and/or alert systemto which the control unit of the drug delivery device 103 is connected.The alarm and/or alert system can be configured to emit a visual and/oraudio alarm and/or alert. The alarm and/or alert system can furtherinclude pre-programmed alarm and/or alert limitations so that when auser modifies any aspect of the system to be outside of the limitations,or the machine itself detects any aspects of the system to be outside ofthe limitations, the module emits an alarm and/or alert.

Turning back to modeling the life cycle of red blood cells, anerythrocyte (i.e., a red blood cell (RBC)) that reaches an age of 120days, has traveled approximately 480 kilometers, making about 170,000circuits, in each cycle enduring osmotic swelling and shrinkage and anumber of deformations while passing through capillaries. Theaccumulated damage to the membrane of the RBC is believed to lead to thedestruction of the cell.

About two decades ago it was observed that erythrocytes actually undergosuicidal death. The suicidal death of erythrocytes is called eryptosis.This term was chosen in order to point out the differences from andsimilarities to apoptosis, the programmed cell death of nucleated cells.Although there exist a number of studies concerning eryptosis, and someparts are understood in great detail, the precise process and mechanismsremain elusive.

Stimulation of eryptosis is governed by a complex signaling network andinvolves activation of cation channels. (F. Lang, C. Birka, S. Myssina,K. S. Lang, P. A. Lang, V. Tanneur, C. Duranton, T. Wieder, and S. M.Huber. Advances in Experimental Medicine and Biology, 559:211-217,2004). Shortly before they are phagocytosed by macrophages, which takesplace mainly in the spleen, one can observe distinctive cell shrinkage,plasma membrane microvesiculation and an exposure of phosphatidylserine(PS) on the cell surface. These PS-exposing erythrocytes are recognized,engulfed and degraded by macrophages. (D. Bratosin, J. Estaquier, F.Petit, D. Arnoult, B. Quatannens, J. P. Tissier, C. Slomianny, C.Sartiaux, C. Alonso, J. J. Huart, J. Montreuil, and J. C Ameisen. CellDeath Differ, 8(12):1143-1156, December 2001).

The normal lifespan and senescence can be either accelerated or delayedby environmental signals. Triggers of eryptosis include osmotic shock,oxidative stress, prostaglandin E2, energy depletion, chlorpromazine,aluminum, mercury, etc. Diseases which are associated with an earlyexpression of PS and thus accelerated eryptosis are, for instance, irondeficiency, sepsis, hemolytic uremic syndrome and sickle-cell anemia.See R. F. Zwaal, P. Comfurius, and E. B. Bevers Cellular and molecularlife sciences, 62:971-988, 2005. On the other hand eryptosis may beinhibited by erythropoietin, adenosine, catecholamines, nitric oxide andactivation of G-kinase. See M. Foller, S. M. Huber, and F. Lang. IUBMBLife, 60:661-668, 2008. Thus, EPO seems not only to increase the numberof circulating erythrocytes by preventing apoptosis of progenitor cellsbut also by prolonging the lifespan of circulating RBC by inhibiting thecation channels. See F. Lang, K. S. Lang, P. A. Lang, S. M. Huber, andT. Wieder Antioxidants & Redox Signaling, 8:1183-1192, 2006; A. B.Schwartz, B. Kelch, L. Terzian, J. Prior, K. E. Kim, E. Pequinot, and B.Kahn. ASAIO transactions, 36:M691M696, 1990. The cation channels arevolume-sensitive, and, after activation of the channel,phosphatidylserine asymmetry is reversed and phosphatidylserine isexposed to the outside of the erythrocyte. Binding of erythropoietin toRBC and inhibition of the channels may contribute to an increase in theerythrocyte number during erythropoietin therapy. See S. Myssina, S. M.Huber, C. Birka, P. A. Lang, K. S. Lang, B. Friedrich, T. Risler, T.Wieder, and F. Lang. Journal of American Society of Nephrology,14:2750-2757, 2003. Surprisingly, it was recently observed thaterythrocytes from erythropoietin overexpressing mice die faster ex vivo.See M. Foeller, R. S. Kasinathan, S. Koka, S. M. Huber, B. Schuler, J.Vogel, M. Gassmann, and F. Lang. American Journal ofPhysiology-Regulatory, Integrative and Comparative Physiology,293:R1127R1134, 2007. This observation leaves room for discussion. Ithas been suggested that: “Possibly, erythropoietin stimulates theexpression of genes in progenitor cells, which foster eryptosis and thusaccelerate erythrocyte death as soon as the erythropoietinconcentrations decline.” See F. Lang, E. Gulbins, H. Lerche, S. M.Huber, D. S. Kempe, and M. Faller. Eryptosis, Cellular Physiology andBiochemistry, 22:373-380, 2008.

Erythropoiesis and its Regulatory Mechanisms

The daily production rate in a healthy adult totals about 200×10⁹ redblood cells under normal conditions. This number can vary due to varyinginternal and environmental conditions. The body has to be able to adaptthe number of erythrocytes to situations of anemia and/or hypoxia (e.g.after bleeding, due to a shortened RBC lifespan, due to an enhancederyptosis, . . . ) as well as to situations of excessive red blood cells(e.g., high altitude dwellers descending to sea level).

The production of new erythrocytes takes place in the bone marrow.Undifferentiated stem cells in the bone marrow commit to the erythroidlineage and undergo a series of proliferations and differentiations. Theearliest stage of erythroid committed cells are called BFU-Es(Burst-Forming Unit Erythroid). Within a few days the BFU-Es mature toCFU-Es (Colony-Forming Unit Erythroid), then they undergo differentstages of erythroblasts, and finally they become reticulocytes. Thereticulocytes are released to the blood stream and within 1-2 days theyappear as mature erythrocytes. See Williams Hematology. Throughout thedescription below, the term progenitor cells is used to refer to BFU-Eand CFU-E cells, and the term precursor cells is used to subsumeproerythroblasts, basophilic erythroblasts, orthochromatophilicerythroblasts and bone marrow reticulocytes.

The primary control of erythropoiesis is governed by the hormoneerythropoietin. EPO is produced mainly by peritubular cells in the renalcortex based on a negative feedback mechanism. The kidneys detect thepartial pressure of oxygen in blood and react by releasing more EPO ifthe oxygen content is too low and vice versa. The concentration of EPOaffects the number of circulating RBC by determining the number of cellsthat mature into erythrocytes. Erythropoietin plays a role in therecruitment of stem cells to the erythroid lineage, prevents apoptosisof progenitor cells and affects the maturation velocity of progenitorand precursor cells. Thus, disturbances in oxygen delivery can beadjusted for by an adaptive resetting of the rate of erythropoiesis.Additionally, there exists a physiological process wherein macrophagesstart to phagocytose young erythrocytes (neocytes), which is calledneocytolysis. Neocytolysis seems to be triggered by a drop in the EPOlevel and helps to regulate situations of excessive red cell mass. SeeRice 2001. Further, recent studies suggest that the concentration of EPOin blood influences the clearance of senescent RBCs and that EPO canprolong the lifespan of erythrocytes by inhibiting cation channels.

Progenitor and Precursor Cells

When a stem cell is committed to the erythroid lineage, it develops intoan erythrocyte. The number of stem cells which enter the differenthematopoietic lineages is determined by interleukines and growthfactors. The immature erythrocyte undergoes a number of changes instructure, appearance and its requirements during the process ofdifferentiation and proliferation. BFU-E cells, the first committederythroid cells, express only a very small number of EPO-receptors(EpoR) and therefore they are almost EPO-independent for their survival.Within a few days they develop into CFU-E cells.

As cells mature to CFU-Es, the number of EpoR on the cell surfaceincreases distinctly, and the cells become absolutely dependent onerythropoietin to prevent them from apoptosis. See H. Wu, X. Liu, R.Jaenisch, and H. F. Lodish. Cell, 83(1):59-67, 1995 (hereinafter “Wu etal.”). Under normal conditions, large numbers of generated CFU-E are notsurviving. After that very EPO-sensitive phase, the density of EpoRsdeclines sharply on early erythroblasts and EpoRs almost disappear atthe stage of orthochromatophilic erythroblasts. See J. P. Greer, J.Foerster, G. M. Rodgers, F. Paraskevas, B. Glader, D. A. Arber, and R.T. Jr. Means. Wintrobe's Clinical Hematology, volume 1. LippincottWilliams & Wilkins, 12th edition, 2009 (hereinafter “Wintrobe's ClinicalHematology”). Under high levels of EPO, the marrow transit time ofprecursor cells is shortened, and high EPO levels result in release ofstill immature reticulocytes, which are referred to as stressreticulocytes. (See Williams Hematology).

One significant difference between progenitor and precursor cells is thesynthesis of hemoglobin. Unlike progenitor cells, all types of precursorcells synthesize heme and the increased demand for iron for this processis reflected by a sharp increase in the number of transferrin receptors(TfR) in proerythroblasts (transferrin receptors are the path of ironuptake). TfRs reach their peak in intermediate erythroblasts and declineafterward and mature bone marrow reticulocytes express only a few TfRs(in comparison).

Neocytolysis

The body has to be able to deal with situations of too few circulatingerythrocytes as well as with excessive red blood cells. Whereas it isvery effective to change the rate of erythropoiesis to adapt the systemto situations when red cell numbers are low, this is not an adequateregulatory mechanism in situations of red cell mass excess. This isbecause a decrease in the rate of erythropoiesis takes relatively longto have a noticeable effect on the number of circulating red bloodcells, due to the long lifespan of erythrocytes.

A mere 20 years ago, preventing apoptosis of red cell progenitors hadbeen thought to be the sole regulating effect of EPO. Because this wouldonly allow for a slow adaption in situations of too many RBC, thecontrol of the body of red cell mass would be coarse. A decline inerythropoietin level results in more progenitor cells dying, but asuppression of the hormone has no effect on maturation of erythroidprecursors. Thus, no decrement in the erythrocyte production would beobservable within one week after EPO has declined. However, studies doneon residents at very high altitude who rapidly descended to sea level,showed a decrease in red cell mass of 10-18% in the first 7-10 days.(See Rice 2001). Persons living in hypoxic environments, like highaltitude, are normally polycythemic. Thus, their hematocrit is too highunder normal conditions and as a consequence the release of EPO issuppressed.

Further investigations of situations where EPO levels are lower thannormal (polycythemic high altitude dwellers, anemia of renal failure,human model based on EPO administration then withdrawal), suggest that asuppression of EPO leads to selective hemolysis of young red bloodcells. (See Rice 2005). This process is called neocytolysis, to stressthe fact that young erythrocytes (i.e., neocytes) are uniquelysusceptible. See C.-C. Chang, Y. Chen, K. Modi, O. G. Awar, C. P.Alfrey, and L. Rice. Journal of Investigative Medicine, 57:650-654,2009; Rice 2005; M. M. Udden, Pickett M. H. Driscoll, T. B., C. S.Leach-Huntoon, and C. P. Alfrey. The Journal of Laboratory and ClinicalMedicine, 125:442-449, 1995.

Neocytolysis is initiated by a fall in erythropoietin levels and in Rice2005 it was shown that, for instance, low doses of EPO, administered tohigh altitude dwellers on descent, completely abrogated neocytolysis. Atthe moment it remains unclear whether it is the rate of decline in EPOlevel or the drop of EPO beneath a certain threshold that acts as atrigger for neocytolysis. See C. P. Alfrey and S. Fishbane. NephronClinical Practice, 106:149-156, 2007 (hereinafter “Alfrey et al.”); Rice1999. The current approach to treatment of anemia with ESAs deviatessignificantly from the normal internal systemic process. Neocytolysis,for instance, contributes to renal anemia. It is precipitated by thepathologic endogenous erythropoietin deficiency of renal disease and theshort bursts in EPO concentration followed by sharp declines in serumlevels due to administration of the hormone (especially duringintravenous (i.v.) administration). It would be desirable to choosedosing schedules such that neocytolysis is minimized or totallyprevented. (See Alfrey et al.; Rice 1999). A mathematical model could bevery useful to define better ESA administration schemes which stimulateerythropoiesis in a more natural way and abrogate neocytolysis. Betterand/or alternative administration schemes can include daily ESAadministration, continuous ESA administration (e.g., via a pump), ormonthly ESA administration.

Iron and Erythropoiesis

Erythropoiesis is a very complex process and even though erythropoietinis the key regulator, there are other proteins (e.g., interleukins,etc.) and substances (e.g., iron, folic acid, vitamin B.sub.12, etc.)that are needed for an optimal erythropoiesis. Iron is a very criticalfactor for an effective production of red blood cells. The metal isindispensable for hemoglobin synthesis and the hemoglobin protein is theactual oxygen transporter in erythrocytes. The molecule makes up about97% of the dry weight of red blood cells. If the body is not able toprovide sufficient iron for the production of erythrocytes, thenimpaired erythropoiesis will result. (S. Finch, D. Haskins, and C. A.Finch. The Journal of Clinical Investigation, 29:1078 1086, 1950;Williams Hematology). In general, a healthy adult has difficultyproviding sufficient iron to support production rates greater than threetimes basal. Higher rates may be possible when administering EPO andiron and in some diseases, e.g. hemochromatosis. (L. T. Goodnough.Nephrology Dialysis Transplantation, 17:14-18, 2002 (hereinafter“Goodnough 2002”); C. A. Finch. Blood. The Journal of American Societyof Hematology, 60(6):1241-1246, 1982) (hereinafter “Finch 1982”). Anundersupply of iron additionally leads to an increase in the number ofhypochromic RBC. Hypochromic cells (i.e., cells in which the amount ofhemoglobin is lower than the normal 26 pg) are small, have a reducedoxygen carrying capacity and are relatively fragile, and thus are likelyto die earlier than normochromic erythrocytes. See D. M. Wrighting andN. C. Andrews. Iron Homeostasis and erythropoiesis. Current Topics inDevelopmental Biology, 82:141-159, 2008 (hereinafter “Wrighting”); A.Loria, L. Sanchez-Medal, R. Lisker, E. De Rodriguez, and J. Labardini.Br J Haematol, 13(3):294-302, May 1967.

Erythroid precursor cells are the most avid consumers of iron in thebody. See P. Ponka and C. N. Lok. The International Journal ofBiochemistry & Cell Biology, 31:1111-1137, 1999 (hereinafter “Ponka1999”). The immature erythrocyte has only a few days to synthesize allthe hemoglobin which the mature cell contains. Hemoglobin is ametalloprotein. The name refers to its special structure. A hemoglobinmolecule consists of four subunits each composed of a globular proteinembedding a heme group, and every heme group in turn contains one ironatom. Erythroid cells rely completely on transferrin receptors to takeup iron. See R. Y. Chan, C. Seiser, H. M. Schulman, L. C. Kuhn, and P.Ponka. European Journal of Biochemistry, 220:683 692, 1994. Unlikeprogenitor cells, all types of precursor cells synthesize heme and theincreased demand for iron for this process is reflected by a sharpincrease in the number of TfR in proerythroblasts. Transferrin receptorsreach their peak in intermediate erythroblasts followed by a decreasewith further maturation.

In precursor cells, heme is essential for maintaining a “normal” numberof TfR. Studies showed that inhibition of heme synthesis stronglyinhibited TfR expression. It seems that there exists a positive feedbackmechanism in which heme promotes a high rate of transferrin receptorsynthesis. A high number of transferrin receptors enhances iron uptakeand that in turn keeps hemoglobin synthesis at high levels. See P.Ponka. Blood, 89:1-25, 1997. There is evidence that precursor cellswhich have a low hemoglobin content are still able to proliferate but donot differentiate and undergo apoptosis. See J. A. Schmidt, J. Marshall,M. J. Hayman, P. Ponka, and H. Beug. Cell, 46:41-51, 1986 (hereinafter“Schmidt”).

1. Iron Homeostasis

Iron plays a vital role in metabolic processes in all living organisms,including oxygen transport, DNA synthesis, and electron transport. It isthe most important essential trace element and amounts to about 35-45and 50-55 mg/kg body weight in adult women and men, respectively. SeeWilliams Hematology; Ponka 1999. For a very good overview of ironhomeostasis and metabolism, see Crichton and Lieu. Homeostasis is thestate of equilibrium (balance between opposing pressures) in the bodywith respect to various functions and to the chemical compositions ofthe fluids and tissues. See Steadman's Medical Dictionary, 26.sup.thEdition, Williams & Wilkins, 1995.

Iron is a key component of many cellular enzymes (e.g. oxidases,catalases, peroxidases, cytochromes, etc.). These enzymes are criticalin many basic cellular processes, including cell proliferation anddifferentiation, DNA and RNA synthesis and electron transport. Further,iron, bound in heme, is essential for oxygen transport. Under normalconditions, total body iron burden is very tightly regulated, becauseexcessive iron leads to tissue damage. Once an atom of iron enters thebody, it is conserved with remarkable efficiency and can remain in thebody for more than ten years. Iron is only lost via loss of cells(epithelial cells, bleeding) and very small amounts are lost via urineand sweat. Unlike other mammals, humans are not able to actively excreteiron. Thus, because of the toxic nature of iron, it is very important tokeep absorption and needs well balanced.

Iron is absorbed from dietary sources via the duodenum. The uptake isregulated by the hormone hepcidin. Under normal conditions, the dailyabsorption of iron amounts to 1-2 mg. In general, in the organism ironis bound to proteins, because free iron results in formation of freeradicals and destroys cells. The iron exporter ferroportin is needed toegress the metal from cells (e.g. enterocytes, macrophages,hepatocytes). Ferroportin is not only the effector of cellular ironexport but also the receptor of the iron regulating hormone hepcidin.See Williams Hematology.

Much of the iron in the human body can be found in circulatingerythrocytes incorporated in hemoglobin (1 ml of packed cells containabout 1 mg of iron). However, most of the iron used for the productionof new erythrocytes comes from hemoglobin recycling and not fromabsorption. When a senescent red cell is phagocytosed, the macrophageinternalizes hemoglobin, degrades it, and exports it back out into bloodcirculation.

Iron is stored in the form of ferritin or hemosiderin, where the latteris more a kind of a long-term storage, because iron within deposits ofhemosiderin is very poorly available. Storage sites can be found mainlyin the liver, the spleen, and the bone marrow but, for instance, smallamounts of hemosiderin can be found in almost any tissue. There aredifferences between men and women in the amount of iron that is stored(men: about 200-400 mg, women: about 1000 mg,). See C. Crepaldi, A.Brendolan, V. Bordoni, M. R. Carta, V. D. Intini, F. Gastaldon, P.Inguaggiato, and C. Ronco. Hemodialysis International, 7(3):216-221,2003 (hereinafter “Crepaldi”).

Iron is carried around in plasma by a transport protein, which is calledtransferrin. Transferrin is produced by the liver. Cells which are inneed for iron express transferrin receptors, and, after transferrinforms a complex with this receptor, iron is transported into the cell.Transferrin exists in three different forms: iron-free (apotransferrin),one iron (monoferric transferrin) and two iron (differic transferrin)atoms bound. The concentration of iron and the amount of transferrinpresent in plasma determine the relative abundance of each form.Although iron bound to transferrin amounts to only about 3 mg of totalbody iron, it is the most important iron pool with the highest rate ofturnover, about 30 mg/day. Approximately 80% of this iron is transportedto the bone marrow and used by erythroid cells. See Lieu. The amount oftransferrin-bound iron is determined by three processes: macrophage ironrecycling, iron storage (hepatic) and duodenal iron absorption. SeeWrighting.

In the last 15 years, understanding of the regulation of systemic ironhomeostasis has changed substantially. The antimicrobial peptidehepcidin (hepatic bactericidal protein), which is secreted by the liver,was identified to be the systemic iron regulator.

Hepcidin inhibits iron transport by binding to the iron exporterferroportin. Hence, it keeps gut enterocytes from secreting iron intocirculation, thereby functionally reducing iron absorption. Further, itprevents iron release from macrophages and blocks the stores as well, byshutting off the means of transport out of cells.

There exist some sensing mechanisms in the liver that gauge the amountof iron needed to support erythropoiesis. Note that hepatocytes in theliver, which secrete hepcidin, also release transferrin into circulationand are one of the major iron storage sites. If the concentration ofiron in stores is too low and/or erythropoiesis is increased, then thelevel of the hormone hepcidin is decreased. Thus, for instance, lowhepcidin concentrations are observed in patients with absoluteiron-deficiency anemia, or anemias with high erythropoietic activity.See L. T. Goodnough, E. Nemeth, and T. Ganz. Blood, 116:4754-4761, 2010(hereinafter “Goodnough 2010”). The consequence of a lower hepcidinlevel is that more iron is taken up via the duodenum and more iron isreleased from macrophages and from the stores. See Crichton; Fleming. Ifiron stores are full and/or erythropoiesis is decreased, then morehepcidin is secreted by hepatocytes. As a result of the binding to theiron exporter ferroportin, dietary iron absorption, macrophage ironrecycling, and release of iron from stores are partly or fully blocked.

Hepcidin is an acute phase-reactant, and, during inflammation, cytokinesstimulate overproduction of hepcidin. Hence, the feedback is alteredduring inflammation, which is a prominent problem in chronic kidneydisease. As a consequence, the body is not able to sufficiently supplydeveloping red blood cells with iron and erythropoiesis is impaired,even though there might be more than enough iron in stores andmacrophages, but it is simply not available. This paradoxical situationis called functional iron deficiency, whereas a situation when the ironcontent of the body is too low (depleted stores) is referred to asabsolute iron deficiency.

2. A Mathematical Model for Erythropoiesis Including Iron Homeostasis

For healthy persons, one can assume that erythropoiesis is alwayssufficiently supplied with iron. This assumption certainly does not holdfor some types of anemia and it is not valid for the majority ofdialysis patients. Therefore, the model developed in this sectionincorporates an extensive iron model in order to expand the possibleapplications to all dialysis patients and to pathologies where anunbalanced iron homeostasis leads to an abnormal erythropoiesis (e.g.,anemia of absolute and functional iron deficiency, hemochromatosis,etc.). The incorporation of an iron submodel and its effects onerythropoiesis has a significant influence on the structure of thepopulation equations.

A reduced reaction of the body to erythropoietin or ESAs is often theresult of an insufficient iron availability in the body. A model whichgathers the balance between developing precursor cells and iron that canbe delivered to the bone marrow for hemoglobin synthesis is needed.Further, iron homeostasis itself is a very complex process. There arecomplex feedback mechanisms which determine how much iron is shiftedaround from one place to another in the body and how much is absorbedvia dietary input. Unfortunately, most of these feedback mechanisms arenot well understood.

A compartmental model for iron is considered and linked to a structuredpopulation model which governs the dynamics of the erythroid cellpopulation. For a detailed explanation of amendments and modificationsin the population equations for the iron model see Section 2.1. InSection 2.2 the derivation of the iron model is explained. For anoverview of how the model developed in this section is organized, seeFIG. 9B.

2.1. Population Model.

The structured population equations are designed to link the iron modelto the cell population model in a physiologically meaningful way. Fivedifferent cell population classes (BFU-E, CFU-E, erythroblasts, marrowreticulocytes, circulating RBCs) are considered. The inclusion of ironmakes it necessary to consider some physiological processes in detail.In many applications of structured population equations, the onlyattribute used to distinguish cells is the cell age. This concept isinsufficient for the population class of the precursor cells, where ironis used to synthesize hemoglobin. Here, cells are characterized by theirhemoglobin content, and the number of tranferrin receptors on the cellmembrane, in addition to the cell age. For mature erythrocytes, cell ageand hemoglobin content are used as attributes.

The two equations describing the progenitor cell populations (BFU-Es andCFU-Es) are solely structured using cell age. These types of cells donot synthesize hemoglobin and, therefore, they are not among the majorconsumers of iron. Hence, these two classes can be described bypopulation equations of the form

${{{\frac{\partial}{\partial t}{u\left( {t,\mu} \right)}} + {{v\left( {E(t)} \right)}{\frac{\partial}{\partial\mu}{u\left( {t,\mu} \right)}}}} = {\left( {{\beta\left( {{E(t)},\mu} \right)} - {\alpha\left( {{E(t)},\mu} \right)}} \right){u\left( {t,\mu} \right)}}},$

where u(t, μ) is the population density of cells at time t with cell ageμ. The functions β and α describe the rate of division and rate ofapoptosis, respectively, of the cells. Further, v describes thematuration velocity.

Precursor cells are the only erythroid cells that synthesize hemoglobin.In order to introduce these phenomena in the model and to keep track ofthe amount of hemoglobin the cells are synthesizing, two additionalattributes are introduced in addition to cell age:

-   -   the number τ of TfRs per cell    -   the amount h of hemoglobin each cell is carrying

Thus, the general population equation for precursor cells is

${{\frac{\partial}{\partial t}{u\left( {t,\mu,h,\tau} \right)}} + {\frac{\partial}{\partial\mu}\left( {{v_{1}\left( {t,\mu,h,\tau} \right)}{u\left( {t,\mu,h,\tau} \right)}} \right)} + {\frac{\partial}{\partial h}\left( {{v_{2}\left( {t,\mu,h,\tau} \right)}{u\left( {t,\mu,h,\tau} \right)}} \right)} + {\frac{\partial}{\partial\tau}\left( {{v_{3}\left( {t,\mu,h,\tau} \right)}{u\left( {t,\mu,h,\tau} \right)}} \right)}} = {{8\beta\left( {t,\mu,{2h},{2\tau}} \right){u\left( {t,\mu,{2h},{2\tau}} \right)}} - \left( {{\beta\left( {t,\mu,h,\tau} \right)} + {\sigma\left( {t,\mu,h,\sigma} \right){{u\left( {t,\mu,h,\tau} \right)}.}}} \right.}$

Here u(t, μ, h, τ) denotes the population density of the cell populationat time t with maturity μ, hemoglobin content h and number oftransferrin receptors τ Further, the functions β and α describe the rateof division and rate of apoptosis of the cells. The functions v₁, v₂, v₃define the rate with which the structural variables μ, h, and τ, arechanging, respectively. For a derivation of the partial differentialequation (1), see below.

Circulating red blood cells carry hemoglobin but they do not expressTfRs for synthesis of the hemoglobin molecule and, therefore, adistinction is made between cells using only cell μ and hemoglobincontent h as structural variables. Hence, the general populationequation reads as follows

${{\frac{\partial}{\partial t}{u\left( {t,\mu,h} \right)}} + {\frac{\partial}{\partial\mu}\left( {{v_{1}\left( {t,\mu,h,} \right)}{u\left( {t,\mu,h} \right)}} \right)} + {\frac{\partial}{\partial h}\left( {{v_{2}\left( {t,\mu,h,} \right)}{u\left( {t,\mu,h} \right)}} \right)}} = {{- {\alpha\left( {t,\mu,h} \right)}}{{u\left( {t,\mu,h} \right)}.}}$

2.1.1. Progenitor Cells: BFU-E and CFU-E Cells. Assumptions:

-   -   1. The number of stem cells, which commit to the erythroid        lineage depends on EPO.    -   2. Cells normally stay in this stage for 8 days (3 days BFU-E        and 5 days CFU-E).    -   3. The number of transferrin receptors on BFU-E and CFU-E cells        is negligible.    -   4. EPO has no effect on the number of divisions or the rate of        apoptosis of BFU-E.    -   5. The proliferation rate of CFU-E cells is constant.    -   6. The rate of apoptosis depends strongly on EPO levels.    -   7. The maturation velocities of BFU-E and CFU-E cells are        constant.

Different kinds of growth factors affect the number of stem cells and invitro studies suggest that high levels of EPO can result in up to 20%more stem cells committing to the erythroid lineage. BFU-E and CFU-Ecells express only a very small number of transferrin receptors on thecell membrane and the influence of iron on these cells is neglected.This is not entirely true, because, like all proliferating cells, theyneed a small amount of iron for division. These amounts are so smallcompared to the consumption of iron for hemoglobin synthesis that theyare neglected. Even in situations of iron deficiency, the assumption isthat there is enough iron for cell divisions of progenitor cellsavailable and only hemoglobin synthesis is impaired. Thus, the followingpopulation equation is obtained for BFU-E cells:

$\left\{ \begin{matrix}{{{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},} \\{{{p\left( {t,0} \right)} = {S_{0}\left( {E(t)} \right)}},} \\{{{p\left( {0,\mu^{p}} \right)} = {p_{0}\left( \mu^{p} \right)}},}\end{matrix} \right.$

where p(t, μ^(p)) is the population density of the cell class at time twith maturity μ^(p), 0≤μ^(p)≤μ_(max) ^(p)=3, t>0. Further, β^(p) is theconstant proliferation rate (Assumption 4) and P₀(μ^(p)) is thepopulation density at t=0. The number S₀(E(t)) of cells committing tothe erythroid lineage is given by a sigmoidal function,

${{S_{0}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + e^{{{- k_{1}}{E(t)}} + c_{1}}} + b_{1}}},$

where E(t) is the EPO concentration at time t and a₁, b₁, c₁ and k₁ arepositive constants and a₁>b₁. The function S₀ increases monotonicallywith increasing EPO concentration (Assumption 1).

The equation for CFU-E cells is:

$\left\{ \begin{matrix}{{{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},} \\{{{q\left( {t,\mu_{\min}^{q}} \right)} = {p\left( {t,\mu_{\max}^{p}} \right)}},} \\{{{q\left( {0,\mu^{q}} \right)} = {q_{0}\left( \mu^{q} \right)}},}\end{matrix} \right.$

where q (t, μ^(q)) is the population density of the CFU-E class at timet with maturity μ^(q), t>0 and 3=μ^(q) _(min)≤μ^(q) _(min)=8. Theconstant β^(q) describes the division rate (Assumption 4) and α^(q)(E(t)) the apoptosis rate which depends on the EPO-concentration(Assumption 5). Moreover, the number of cells leaving the BFU-E cellstage and entering the CFU-E cell stage per unit time requires theboundary condition q(t, μ^(q) _(min))=p(t, μ^(q) _(max)). Finally, q₀(μ^(q)) is the population density at t=0.

Again, a sigmoidal function is used to describe the dependence of therate of apoptosis on EPO,

${{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + e^{{k_{2}{E(t)}} - c_{2}}} + b_{2}}},$

where E(t) is the EPO concentration and a₂, b₂, c₂ and k₂ are positiveconstants, a₂>b₂. The function α^(q) is per definition a monotonicallydecreasing function. Thus, a higher level of EPO causes more cells tosurvive.

Finally, because of Assumption 7, it is assumed that cell age equals theactual age of the cells, i.e., v^(p)=v^(q)≡1, ∀≥0.

2.1.2. Precursor Cells: Erythroblasts and Marrow Reticulocytes.

Assumptions:

-   -   8. The class erythroblasts consists of all cell stages from        proerythroblast to orthochromatophilic erythroblast.    -   9. Cells normally stay in this stage for 6-8 days (5 days        erythroblasts and 1-3 days marrow reticulocytes).    -   10. Under high levels of EPO, so called “stress reticulocytes”        are released.    -   11. All types of precursor cells synthesize hemoglobin.    -   12. The rate at which the amount of hemoglobin in precursor        cells changes, depends on the number of TfRs on the cell surface        and the currently available amount of iron.    -   13. The rate of change of TfRs in precursor cells depends on the        cell age and the hemoglobin content of the cell.    -   14. Cells entering the erythroblast class may be equipped with        different numbers of TfRs, i.e., TfRs need not be uniformly        distributed among the erythroblast population at time t=0.    -   15. Hemoglobin is conserved during cell division and a mother        cell divides into 2 daughter cells carrying equal amounts of        hemoglobin.    -   16. TfRs are conserved during division and a mother cell divides        into 2 daughter cells with equal numbers of TfRs.    -   17. The maturation velocity of erythroblasts is constant.    -   18. The proliferation rate of erythroblasts is constant.    -   19. The rate of apoptosis of erythroblasts depends on cell age        and hemoglobin.    -   20. The maturation velocity of reticulocytes depends on EPO.    -   21. The rate of apoptosis of reticulocytes depends on cell age        and on hemoglobin.    -   22. Reticulocytes mature but do not proliferate.    -   23. Iron that enters a precursor cell is incorporated into        hemoglobin, i.e., no free iron in the cell is considered.

The number of TfRs increases sharply in the proerythroblasts and reachesits peak in intermediate erythroblasts before declining again(proerythroblasts: 300,000 TfRs/cell; basophilic erythroblasts: 800,000TfRs/cell; mature reticulocytes: 100,000 TfRs/cell). See Wintrobe'sClinical Hematology. This reflects the increased demand of the cells foriron to synthesize hemoglobin Almost all iron contained in the cells issynthesized to hemoglobin, leaving mature erythrocytes with only anegligible amount of non-heme iron (Assumption 23). Depending on thematurity of cells different numbers of TfRs are to be expected(Assumption 13). Cells that are saturated with hemoglobin reduce thelevel of receptors. Further, there is evidence that in contrast to othercells, in the precursor cells the number of transferrin receptors isalso controlled by a positive feedback involving the amount of alreadysynthesized hemoglobin in the cell. In vitro hemoglobin was shown to beessential for these cells to maintain a normal number of TfRs. See Ponka1999.

For the erythroblasts population class, the following PDE is obtained

$\left\{ \begin{matrix}{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} +} \\{{v_{2}^{r}\left( {\tau,{a_{1}(t)}} \right){\frac{\partial}{\partial h^{r}}r}\left( {t,\mu^{r},h^{r},\tau^{r}} \right)},} \\{{+ {\frac{\partial}{\partial\tau^{r}}\left( {{v_{3}^{r}\left( {\mu^{r},h^{r},\tau^{r}} \right)}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} \right)}} = {{8\beta{u\left( {t,\mu^{r},{2h^{r}},{2\tau^{r}}} \right)}} -}} \\{{\left( {\beta + {\alpha^{r}\left( {\mu^{r},h^{r}} \right)}} \right){r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}},} \\{{{r\left( {t,\mu_{\min}^{r},0,\tau^{r}} \right)} = {{q\left( {t,\mu_{\max}} \right)} \cdot {\eta\left( \tau^{r} \right)}}},} \\{{if}\left( {\mu^{r},h^{r},\tau^{r}} \right){is}a{boundary}{point}} \\{{{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)} = 0},} \\{{{{and}h^{r}} > 0},} \\{{{r\left( {0,\mu^{r},h^{r},\tau^{r}} \right)} = {r_{0}\left( {\mu^{r},h^{r},\tau^{r}} \right)}},}\end{matrix} \right.$

where r(t, μ^(r), h^(r), τ^(r)) is the population density of theerythroblast class at time t with maturity μ^(r), hemoglobin contenth^(r) and number of TfRs τ^(r), 8=μ_(min) ^(r)≤μ^(r)≤μ_(max) ^(r)=13,0≤h^(r)≤h_(max) ^(r), τ_(min) ^(r)≤τ^(r)≤τ_(max) ^(r), t>0. Thepopulation density at t=0 is given as (μ^(r), h^(r), τ^(r)). Assumption17 allows setting v_(μ) ^(r)≡1, i.e., the cell age of the cell coincideswith the actual age of the cell. The contact rate between the TfRs of acell and transferrin molecules in plasma carrying iron is governed bythe law of mass action:

d ₁τ^(r)(t)a ₁(t),

where d₁ is a rate constant, τ^(r)(t) is the number of TfRs on the cellmembrane at time t and a₁(t) is the concentration of iron in plasma attime t. Thus, the amount of hemoglobin in a cell evolves according to

v ₂ ^(r)(τ^(r) ,a ₁)=d ₂τ^(r) a ₁,  (2)

where d₂ is a rate constant that is different from d₁ (note that fouratoms of iron are needed for one hemoglobin molecule). The change of thenumber of TfRs on the cell surface is described by a function v_(τ)^(r), depending on the cell age μ^(r) and the amount of hemoglobin h^(r)(Assumption 13). The function v_(τ) ^(r) is chosen such that itincreases in the beginning and then decreases when the cells are “fullysaturated” with hemoglobin.

${{v_{3}^{r}\left( {\mu^{r},h^{r},\tau^{r}} \right)} = {{a_{3}{\tau^{r}\left( {h_{*}^{r} - {b_{3}\left( h^{r} \right)}^{2}} \right)}} - \frac{c_{3}{\max\left( {\tau^{r},0} \right)}\left( {\mu^{r} - \mu_{\max}^{q}} \right)}{1 + h^{r}}}},$

where a₃, h*^(r), b₃ and c₃ are constants and μ^(q) _(max) is themaximal maturity of progenitor cells. The second term assures that cellswhich progress in maturation, but do not accumulate sufficient iron,lose transferrin receptors. The term max(τ^(r), 0) is needed toguarantee that the number of TfRs can not become negative. The maximumfunction is not differentiable, therefore it is approximated by thefollowing function

smax(τ^(r))=0.5τ^(r)+0.5√{square root over ((τ^(r))²+0.01)}.

The rate of division β^(r) is constant (Assumption 15). Note that theterm 8βu(t, μ^(r), 2h^(r), 2τ^(r)) in the equation arises because ofAssumptions 15 and 16. For further details see below. To model the rateof apoptosis α^(r), a decrease in expression of TfRs is taken intoaccount and leads to a decrease in iron uptake and, thus, lesshemoglobin is synthesized. Cells containing small amounts of hemoglobinarrest in differentiation (proliferation is not affected) and this leadsto premature cell death. See Schmidt Therefore, it is assumed that α^(r)depends on μ^(r) and h^(r) (Assumption 19)

$\begin{matrix}{{{\alpha^{r}\left( {\mu^{r},h^{r}} \right)} = \frac{{b_{4}\left( {h^{**} - h^{r}} \right)}^{2}}{1 + \left( {\mu^{**} - \left( {\mu^{r} - \mu_{\max}^{q}} \right)} \right)^{2}}},} & (3)\end{matrix}$

where b₄, h**, μ**>0 are constants. The parameters h** and μ** can, forinstance, be chosen as

h**=h _(max) ^(r),μ**=μ_(max) ^(r)−μ_(max) ^(q),  (4)

where h^(r) _(max) is the maximal hemoglobin amount of precursor cellsand μ^(r) _(max) is the maximal maturity of precursor cells.

In principle, the boundary conditions arise as a result of the fact thatcells from the precedent class leave the class when they reach themaximum maturity and enter the next population class. Hence, there is aconstant flux of cells from the CFU-E class to the erythroblasts class.However, the passage is more complicated because there is a flux from apopulation class with only one structural variable to a population classwith three attributes. Thus, it has to be carefully considered at whichpoints of the boundary the cells are entering. It is useful to keep thephysiological situation in mind. Progenitor cells do not synthesizehemoglobin. Hence, there are no cells entering the erythroblastspopulation class with h>0. On the other hand, it is not natural to applythe same assumption for transferrin receptors. For although the numberof TfRs on progenitor cells is neglected, in fact, very mature CFU-Eexpress a number of TfRs. Further, it seems reasonable to assume that atthe time when they become proerythroblasts not all cells incorporateexactly the same amount of TfRs on the surface, but that there is somedistribution of this attribute among the cells (Assumption 14). Thus,the concept for the boundary includes, that

∫_(τ_(min))^(τ_(max))η(τ)dτ = 1,

where, for instance, η(τ) can be chosen such that the structuralvariable t is normally distributed with a certain mean value andvariance.

The revisited population equation for marrow reticulocytes reads

$\left\{ \begin{matrix}{{\frac{\partial}{\partial t}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} + {{v_{1}^{s}\left( {E(t)} \right)}{\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}} +} \\{{{v_{2}^{s}\left( {\tau^{s},{a_{1}(t)}} \right)}{\frac{\partial}{\partial h^{s}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}},} \\{{{+ {\frac{\partial}{\partial\text{?}}\left( {{v_{3}^{s}\left( {\mu^{s},h^{s},\tau^{s}} \right)}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} \right)}} = {{- {\alpha^{s}\left( {\mu^{s},h^{s}} \right)}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}},} \\{{{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\min}^{s},h^{s},\tau^{s}} \right)}} = {r\left( {t,\mu_{\max}^{r},h^{r},\tau^{r}} \right)}},} \\{{{s\left( {0,\mu^{s},h^{s},\tau^{s}} \right)} = {s_{0}\left( {\mu^{s},h^{s},\tau^{s}} \right)}},}\end{matrix} \right.$ ?indicates text missing or illegible when filed

where s(t, μ^(s), h^(s), τ^(s)) is the population density of the marrowreticulocytes class at time t with maturity μ^(s), amount of hemoglobinh^(s) and number of TfRs τ^(s),

13=μ_(min) ^(s)≤μ^(s)μ_(max) ^(s)=15, h_(min) ^(s)≤h^(s)≤h_(max) ^(s),τ_(min) ^(s)≤τ^(s)≤Σ_(min) ^(s), t>0. The boundary conditionv^(s)(E(t))s(t, μ_(min) ^(s), h^(s), τ^(s))=r(t, μ_(max) ^(r), h^(r),τ^(r)) describes the flux of cells leaving the erythroblast class andentering the reticulocytes population class carrying a certain amount ofhemoglobin and expressing a certain number of TfRs. Finally, s₀(μ^(s),h^(s), τ^(s)) is the population density at t=0.

A sigmoid function is used to describe the changes in the maturationvelocity v_(μ) ^(s),

${{v_{1}^{s}\left( {E(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + e^{{{- k_{5}}{E(t)}} + c_{5}}} + b_{5}}},$

where E(t) is the EPO concentration at time t and a₅, b₅, c₅ and k₅ aresome positive constants with a₅>b₅. The maturation velocity increaseswith a rising concentration of EPO (Assumption 10). Note, a slowermaturation velocity causes the cells to reach the maximum cell age μ^(s)_(max) at a later point, whereas a faster maturation velocity shortensthe transit time, i.e., the cells reach μ^(s) _(max) earlier. The changeof the amount of hemoglobin v_(h) ^(s) in a cell, the development ofTfRs v_(τ) ^(s) on the cell surface and the rate of apoptosis α_(s) aredescribed with the same functions as for the erythroblast class, i.e.

v ₂ ^(s)(τ,a ₁)=v ₂ ^(r)(τ,a ₁),v ₃ ^(s)(μ,h,τ)=v ₃^(r)(μ,h,τ),α^(s)(μ,h)=α^(r)(μ,h).

Further, β^(s)≡0 because of Assumption 22.

So far, Assumption 19 has not been accounted for. The assumption can beapplied in the following way: If E(t)>E* the reticulocyte class isskipped and erythroblasts with maximal maturity directly enter theerythrocytes class.

2.1.3. Erythrocytes.

Assumptions:

-   -   24. Erythrocytes and blood reticulocytes are subsumed in one        class.    -   25. Cells mature but do not proliferate.    -   26. There is a fixed random daily break-down of red blood cells        (not to be confused with loss of erythrocytes due to        senescence).    -   27. A drop in EPO concentration beneath a threshold precipitates        neocytolysis.    -   28. Erythrocytes with age between 14-21 days are likely to be        affected by neocytolysis.    -   29. Erythrocytes lose hemoglobin during senescence.    -   30. The expected life span of normochromic cells in healthy        persons is 120 days, but can decrease significantly in uremic        patients.    -   31. Cells in which hemoglobin content drops beneath a certain        threshold or which reach the maximum cell age are phagocytosed.

Up to this point, it has been assumed that red blood cells carryapproximately the same amount of hemoglobin. This assumption is notvalid in many pathologies, different types of anemia, and the majorityof dialysis patients, where a wide distribution of corpuscularhemoglobin can be observed. Additionally, it is desirable to account forhemoglobin which is lost during senescence. This loss, which amounts toabout 15-20%, is a non-negligible amount of iron that flows back to themacrophages. See S. C. Gifford, J. Derganc, Y. Shevkoplyas, S. S.Yoshida., and M. W. Bitensky, British Journal of Haematology,135:395-404, 2006; F. L. Willekens, F. H. Bosch, B.Roerdinkholder-Stoelwinder, Y. A. Groenen-Dopp, and J. M. Were, EuropeanJournal of Haematology, 58:246-250, 1997; J. M. Werre, F. Willekens, F.H. Bosch, L. D, de Haans, S. G, van der Vegt, A. G, van den Bos, and G.J. Bosman, Cellular and Molecular Biology 50 (2), 139-145, 2004; F. L.Willekens, B. Roerdinkholder-Stoelwinder, Y. A. Groenen-Dopp, H. J. Bos,G. J. Bosman, A. G. van den Bos, A. J. Verkleij, and J. M. Were, Blood,101:747-751, 2003. Further, the binding of oxygen is pivotal fordelivery of oxygen by hemoglobin, and thus a reduction of 15-20% in thehemoglobin mass will have a noticeable impact on the oxygen carryingcapacity. Hence, there are several good reasons that suggest includinghemoglobin content h^(m) as a structural variable. Further, theattribute h^(m) allows one to relatively easily compute the amount ofiron which is set free when an erythrocyte is phagocytosed. Moreover, itenables tracking of how much iron is set free daily because of sheddingof hemoglobin containing vesicles.

Altogether, the following equation is obtained for erythrocytes:

$\left\{ \begin{matrix}{{{{\frac{\partial}{\partial t}m}\left( {t,\mu^{m},h^{m}} \right)} + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}} + {{v_{2}^{m}\left( \mu^{m} \right)}{\frac{\partial}{\partial h}{m\left( {t,\mu^{m},h^{m}} \right)}}}},} \\{{= {{- {\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},} \\{{{m\left( {t,\mu^{m},h_{\min}^{m}} \right)} = 0},} \\{{m\left( {t,0,h^{m}} \right)} = \left\{ \begin{matrix}{v^{s}\left( {{{\left( {E(t)} \right)\text{?}{s\left( {t,\mu_{\max}^{s},h^{s},\tau^{s}} \right)}d\tau^{s}{for}E} \leq E^{*}},} \right.} \\{{{\text{?}r\left( {t,\mu_{\max}^{r},h^{r},\tau^{r}} \right)d\tau^{r}{for}E} > E^{*}},}\end{matrix} \right.} \\{{{m\left( {0,\mu^{m},h^{m}} \right)} = {m_{0}\left( {\mu^{m},h^{m}} \right)}},}\end{matrix} \right.$ ?indicates text missing or illegible when filed

where the m(t, μ^(m), h^(m)) is the population density for theerythrocyte class at time t with maturity μ^(m) and cell hemoglobinh^(m), 0=μ_(min) ^(m)≤μ^(m)≤μ_(max) ^(m)=120, h_(min) ^(m)≤h^(m)≤h_(max)^(m), and t≥0 and m₀(μ^(m), h^(m)) is the population density at t=0.Moreover, m(t, 0, h^(m)) is defined by the number of reticulocytes orthe number of erythroblasts entering the blood stream carrying a certainamount of hemoglobin h^(m), respectively. To describe the corpuscularhemoglobin decay, a function v_(k) ^(m) is used, of the form

v ₂ ^(m)(μ^(m))=b ₆(μ^(m)−μ*^(m))² −c ₆,

with b₆, c₆, μ*^(m)>0 constant. Choosing appropriate constants b₆, c₆and μ*^(m) allows for a wide range of possible forms of the functionh^(m)(t), with

$\frac{{dh}^{m}(t)}{dt} = {{v_{2}^{m}\left( {\mu^{m}(t)} \right)}.}$

The mortality rate α^(m)(E(t), μ^(m), h^(m)) is the sum of three terms:

-   -   1. The rate of cells dying due to senescence, which is assumed        to happen when the cell reaches the maximal cell age μ^(m)        _(max), or the minimal hemoglobin content h^(m) _(min)        (Assumption 30).    -   2. Cells dying because of random break-down.    -   3. Cells dying because of neocytolysis.

Altogether,

${\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)} = \left\{ \begin{matrix}{{{{\gamma^{m}\left( {\mu^{m},h^{m}} \right)} + {\min\left( {\frac{c_{E}}{{E(t)}^{k_{E}}},b_{E}} \right){for}{E(t)}}} < \text{?}},{\mu_{\min}^{m,n} \leq \mu_{\max}^{m,n}},} \\{\gamma^{m}\left( {\mu^{m},h^{m}} \right){otherwise}}\end{matrix} \right.$ ?indicates text missing or illegible when filed

where b_(E), c_(E), k_(E) are positive constants, τ_(E) is the thresholdbeneath which neocytolysis is triggered (Assumption 27). Further,[μ_(min) ^(m,n),μ_(max) ^(m,n)]=[14,21] is the age-interval during whichcells are effected by neocytolysis (Assumption 28) and t≥0. Since it isassumed that a cell is phagocytosed when it reaches the maximal cell ageμ_(max) ^(m) or the minimal cell hemoglobin h_(min) ^(m), the mortalityrate γ^(m)(μ^(m), h^(m)) needs to be chosen such that γ^(m)(μ^(m),h^(m))=α_(rand) ^(m) for μ_(min) ^(m)≤μ^(m)≤μ_(max) ^(m)−δ, h_(min)^(m)+δ≤h^(m)≤h_(max) ^(m) with δ>0 sufficiently small and ∫_(μ) _(max)_(−δ) ^(μ) ^(max) ∫_(h) _(min) ^(h) ^(min+δ) α^(m)(μ, h)dμdh=∞. Here, aa_(nd) denotes the random breakdown of the cells (Assumption 26). Apossible choice for γ^(m)(μ^(m), h^(m)) is

${\gamma^{m}\left( {\mu^{m},h^{m}} \right)} = \left\{ \begin{matrix}{\alpha_{r}^{m}} & {{{for}\mu^{m}{\epsilon\left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack}},} \\ & {{h^{m} \in \left\lbrack {{h_{\min}^{m} + \delta},h_{\max}^{m}} \right\rbrack},} \\{\gamma_{\mu}^{m}\left( \mu^{m} \right)\gamma_{h}^{m}\left( h^{m} \right)} & {{{for}{\mu^{m}\left( {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)}},} \\ & {{h^{m} \in \left( {h_{\min}^{m},{h_{\min}^{m} + \delta}} \right)},} \\\infty & {{{{for}\mu^{m}} \geq \mu_{\max}^{m}},{h^{m} \leq h_{\max}^{m}},}\end{matrix} \right.$ where${{\gamma_{\mu}^{m}\left( \mu^{m} \right)} = \frac{3\alpha_{rand}^{m}{\delta}^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}}},$and${\gamma_{h}^{m}\left( h^{m} \right)} = {\frac{3{\delta}^{2}}{\left( h^{m} \right)^{2} - {2\left( {h_{\min}^{m} + \delta} \right)h^{m}} + {\left( {h_{\min}^{m} + {2\delta}} \right)h_{\max}^{m}}}.}$

In case of bleeding an additional term

a_(bleed) ^(m)(t) needs to be added. What this term looks like dependson whether the bleeding lasts minutes or hours or is prolonged forseveral days or weeks.

2.1.4. Erythropoietin.

Assumptions:

-   -   32. Release of EPO is controlled by a negative feedback        mechanism according to the oxygen content.    -   33. Oxygen carrying capacity is directly proportional to the        amount of hemoglobin of circulating red blood cells.    -   34. The degradation rate of EPO is constant.    -   35. There is a slight delay in reaction of the EPO production        rate to the number of RBCs but this is negligible compared to        the duration of development of erythrocytes.

The equation which describes the release of EPO by the kidneys is

${E_{in}^{end}(t)} = {\frac{a_{7} - b_{7}}{1 + \text{?}}b_{7,}}$?indicates text missing or illegible when filed

where

${\overset{\sim}{H}{b(t)}} = \frac{{Hb}(t)}{TBV}$

is the hemoglobin concentration. Here, TBV is the total blood volume andHb(t)=∫_(k) _(min) _(m) ^(k) ^(max) ^(m) m(t, μ^(m), h^(m))dh^(m) is thetotal amount of hemoglobin in blood.

Furthermore, a₇, b₇, c₇ and k₇ are positive constants with a₇>b₇.

The dynamics of the endogenous EPO concentration E^(end)(t) aredescribed by the following ordinary differential equation:

${{{\overset{.}{E}}^{end}(t)} = {{\frac{1}{TBV}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},$

where E^(end)(t) is the endogenous EPO concentration, E_(in) ^(end) isthe amount of EPO released by the kidneys and c_(deg) ^(end) describesthe degradation rate of endogenous EPO. The change in the concentrationE^(ex)(t) of an ESA reads

${{{\overset{.}{E}}^{ex}(t)} = {{\frac{1}{TBV}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{E^{ex}(t)}}}},$

where E_(in) ^(ex)(t) is the rate at which the artificial hormone isadministered c_(deg) ^(ex) is the rate with which the exogenous hormoneis degraded. In intravenous administration, the total amount of theagent is injected into a vein, within a very short time interval. Inthis case E_(in) ^(ex)(t) can be approximated by E₀ ^(ex)(t)δ_(t0)(t),where E₀ ^(ex) is the amount of artificial hormone administered andδ_(t0)(t) is the Dirac delta impulse located at t₀, the time when theadministration takes place. In addition, the degradation rate forexogenous EPO c_(deg) ^(ex) differs from the one for endogenous EPO andvaries according to the kind of ESA administered. The overallconcentration of EPO in plasma consists of the naturally producederythropoietin in the body and the administered ESA

E(t)=E ^(ex)(t)+E ^(end)(t).

3. The Iron Submodel.

Assumptions:

-   -   37. Five iron pools are considered within the body: a plasma        pool, iron in precursor cells, iron in erythrocytes, iron in the        macrophages and a storage pool.    -   38. Other iron in the body is ignored.    -   39. Hepcidin regulates iron homeostasis.    -   40. Stressed erythropoiesis and depleted stores decrease the        amount of hepcidin released into plasma.    -   41. Large iron stores and inflammation increase hepcidin levels.    -   42. The amount of iron in plasma is proportional to the amount        of transferrin molecules carrying iron. (Free iron is        neglected.)    -   43. No distinction is made between monoferric and differic        transferrin.    -   44. Iron is lost via loss of cells.    -   45. Iron lost via urine and sweat is neglected.

Iron is a part of hemoglobin which provides the means of O₂ transport tothe tissue. The rate at which iron can be released to plasma fromexisting stores and macrophages may easily limit the rate of irondelivery for hemoglobin synthesis. Therefore, a strong relationshipbetween impaired erythropoiesis and total body iron burden exists. Anexcellent description of iron homeostasis and its effects can be foundin Crichton.

Iron is a substance which is very tightly controlled from the body asiron overload is toxic. Iron homeostasis can be only achieved by thecontrol of absorption because the human body, unlike other mammals, isnot able to influence the excretion of iron.

Under normal circumstances the average daily loss of iron is very small(men: 1 mg, women: 2 mg, because of menstruation or pregnancy).

For the iron submodel, five dynamic iron pools are considered within thebody: a plasma pool, iron in precursor cells, iron in erythrocytes, ironin the macrophages and a storage pool. Hence, ferritin and hemosiderinstores are combined. Other iron in the body is ignored. The amount ofiron lost via urine and sweat is neglected. Thus, in this model,excretion of iron takes place only when cells are lost. On one hand,this can be due to loss of RBCs, i.e., due to external bleeding, and onthe other hand there is a daily loss of epithelial cells, which amountfor about 1 mg/day. Further, it is assumed that only the precursor cellsin the bone marrow uptake iron and ineffective erythropoiesis depends onthe amount of iron which can be provided for erythropoiesis. Moreover,the model accounts for the loss of iron from erythrocytes duringsenescence.

The general form of the compartmental iron model is as follows:

${{\frac{d}{dt}a_{1}} = {{{k_{41}(H)}a_{4}} + {k_{{iv},1}a_{iv}} + k_{{gastro},1} + {(H)a_{gastro}} + {{k_{51}(H)}a_{5}} - {k_{51}a_{1}} - {k_{12}a_{1}}}},$${{\frac{d}{dt}a_{2}} = {{k_{12}a_{1}} - f_{24} - f_{23}}},$${{\frac{d}{dt}a_{3}} = {f_{23} - f_{34} - f_{3,{out}}}},$${{\frac{d}{dt}a_{4}} = {f_{24} - {f_{34}a_{3}} - {{k_{41}(H)}a_{4}} - {k_{45}a_{4}}}},$${\frac{d}{dt}a_{5}} = {{k_{15}a_{1}} + {k_{45}a_{4}} - {{k_{51}(H)}a_{5.}}}$

Note, in the above equations, the dependence of particular components ont was not explicitly stated to simplify notation. Furthermore, a_(i),i=1, . . . , 5, denotes the amount of iron in the compartment i andk_(ij), i, j=1, . . . , 5 are the transfer rates from compartment i tocompartment j, i.e., k_(ij)a_(i) is the rate at which iron is movingfrom compartment i to compartment j. Further, f_(ij) denote the rates atwhich iron is transferred from compartment i into compartment j, when itis not of the simple form k_(ij)a_(i); f_(3,out) is the amount of ironwhich is lost by bleeding. Additionally, a_(iv) denotes the amount ofiron intravenously administered, whereas a_(gastro) is the amount ofiron in the duodenum. Finally, k_(iv,1) and k_(gastro,1)(H) are thecorresponding flow rates. The rate k_(gasto,1) can be described by adecreasing sigmoidal function:

${{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{8} - b_{8}}{1 + e^{{{- k_{8}}{H(t)}} + c_{8}}}b_{8}}},$

where a₈, b₈, c₈ and k₈ are positive constants with a₈>b₈. The functionH(t) is the solution of the differential equation (5).

Since the time needed to administer iron is very short compared to othertime constants in the system, the following can be taken

k _(iv,1) a _(iv) =a _(iv,total)δ_(t) ₀ .

Here a_(iv,total) is the total amount of iron administered and δ_(t0) isthe Dirac delta function located at t₀.

3.1. Hepcidin Feedback.

The regulator of iron homeostasis is hepcidin, a peptide hormoneproduced by the liver. Hepcidin negatively regulates the availability ofiron in plasma by binding to ferroportin on the membranes ofiron-exporting cells. Hepcidin binding induces the endocytosis andproteolysis of ferroportin, and thereby decreases the delivery of ironto the plasma. See Nemeth E., Tuttle M. S., Powelson J., Vaughn M. B.,Donovan A., Ward D. M., Ganz T., and Kaplan J., Science 306 pp.2090-2093 (2004). Ferroportin is the iron exporter required for ironegress from iron exporting cells, as for instance, enterocytes andmacrophages. Hepcidin acts as a medium to withhold iron from certaininvasive bacteria, as they are not able to proliferate under conditionsof insufficient iron supply. The dynamics of the hepcidin concentrationin plasma is calculated according to (5)

${{\frac{d}{dt}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{H(t)}}}},$

where H_(prod)(t) is the rate at which hepcidin, produced by thehepatocytes, is released into the plasma at time t. TBV is the totalblood volume, and c^(H) _(deg) is the rate at which hepcidin isdegraded. The production rate H_(prod)(t) of hepcidin is obtained as theresult of a feedback loop involving the total precursor cell population,the amount of iron stored and the state of inflammation of the patient

H _(prod)(t)=ƒ(U(t),a ₅(t),i(t)).

Here

${{U(t)} = {{\overset{\mu_{\max}^{r}}{\int\limits_{\mu_{\min}^{r}}}{\overset{h_{\max}^{r}}{\int\limits_{h_{\min}^{r}}}{\overset{\tau_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}}{{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\overset{\mu_{\max}^{s}}{\int\limits_{\mu_{\min}^{s}}}{\overset{h_{\max}^{s}}{\int\limits_{h_{\min}^{s}}}{\underset{\tau_{\min}^{s}}{\int\limits^{\tau_{\max}^{s}}}{{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}}},$

where U(t) describes the total precursor cell population(erythroblasts+reticulocytes). Further, a₅(t) is the amount of ironstored and i(t) is the inflammation status at time t. The feedback ismodeled using a monotonically increasing sigmoidal function

${{H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + e^{{{- k_{H}}{C(t)}} + c_{H}}} + b_{H}}},$

where a_(H), b_(H), c_(H) and k_(H) are positive constants witha_(H)>b_(H). The function

C(t)=c _(U) U(t)+c _(a) ₅ a ₅(t)+c _(i) i(t).

is an auxiliary equation that increases if storage amount orinflammation status rises, and decreases when more RBCs are produced.Here, c_(U), c_(a5) and c_(i) are positive constants.

3.2. Detailed Description of the Iron Compartments. 3.2.1. Plasma(Compartment 1).

The rate at which iron enters erythroid cells at time t is given by

${d_{0}{a_{1}(t)}\text{⁠}\left( {{\overset{\mu_{\max}^{r}}{\int\limits_{\mu_{\min}^{r}}}{\overset{h_{\max}^{r}}{\int\limits_{h_{\min}^{r}}}{\overset{\tau_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}}{\tau^{r}r\left( {t,\mu^{r},h^{r}} \right)d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\overset{\mu_{\max}^{s}}{\int\limits_{\mu_{\min}^{s}}}{\overset{h_{\max}^{s}}{\int\limits_{h_{\min}^{s}}}{\underset{\tau_{\min}^{s}}{\int\limits^{\tau_{\max}^{s}}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}} \right)},$

where r(t, μ^(r), τ^(r), h^(r)) and s(t, μ^(s), τ^(s), h^(s)) denote thepopulation densities of hemoglobin synthesizing cells, erythroblasts andreticulocytes, respectively. At the same time this is the rate at whichiron leaves the plasma compartment and enters the precursor cellcompartment.

Only a very small amount of iron can be found in this pool (about 3 mg;≤1%) bound to transferrin, but it is a very important compartmentbecause of its rapid turnover rate. See Finch 1982. Iron normally turnsover at least 10 times each day. See Williams Hematology. Inflow intothe plasma compartment comes from macrophages, the storage pool, ironabsorbed through the duodenum and via intravenously administered iron.Outflow leaves the plasma compartment to the precursor cells and thestorage compartment. The flow rates k₄₁, k_(gastro,1) and k₅₁ aredependent on the hepcidin concentration H(t) at time t. They decreasewhen the hepcidin level rises and vice versa. Inflammation and full ironstores increase the release of hepcidin from the hepatocytes in theliver, whereas a stressed erythropoiesis and depleted stores suppressthe production of hepcidin. The outflow to the precursor cellcompartment k₁₂a₁ is altered by the amount of iron needed from the cellsto synthesize hemoglobin. As erythroid cells are only able to take upiron through transferrin receptors, their need for iron can beassociated with the number of receptors expressed on precursor cells.Thus, the rate k₁₂ is given by

$k_{12} = {d_{0}{\left( {{\overset{\mu_{\max}^{r}}{\int\limits_{\mu_{\min}^{r}}}{\overset{h_{\max}^{r}}{\int\limits_{h_{\min}^{r}}}{\overset{\tau_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}}{\tau^{r}r\left( {t,\mu^{r},h^{r}} \right)d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\overset{\mu_{\max}^{s}}{\int\limits_{\mu_{\min}^{s}}}{\overset{h_{\max}^{s}}{\int\limits_{h_{\min}^{s}}}{\underset{\tau_{\min}^{s}}{\int\limits^{\tau_{\max}^{s}}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}} \right).}}$

For an explanation of the various terms see the statements concerningerythroblasts in Subsection 2.1.2. Hence, the change of amount of ironin the plasma compartment is

${\frac{d}{dt}{a_{1}(t)}} = {{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {k_{{gastro},1}\left( {{{H(t)}{a_{gastro}(t)}} + {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{0}\left( {{\overset{\mu_{\max}^{r}}{\int\limits_{\mu_{\min}^{r}}}{\overset{h_{\max}^{r}}{\int\limits_{h_{\min}^{r}}}{\overset{\tau_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}}{\tau^{r}r\left( {t,\mu^{r},h^{r}} \right)d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\overset{\mu_{\max}^{s}}{\int\limits_{\mu_{\min}^{s}}}{\overset{h_{\max}^{s}}{\int\limits_{h_{\min}^{s}}}{\underset{\tau_{\min}^{s}}{\int\limits^{\tau_{\max}^{s}}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}} \right)}{a_{1}(t)}}} \right.}}$

The transfer rate k₄₁(H(t)) and k₅₁(H(t)) can be described bymonotonically decreasing sigmoidal functions

$\begin{matrix}{{{k_{41}\left( {H(t)} \right)} = {\frac{a_{9} - b_{9}}{1 + e^{{{- k_{9}}{H(t)}} + c_{9}}} + b_{9}}},} & (6)\end{matrix}$ and $\begin{matrix}{{{k_{51}\left( {H(t)} \right)} = {\frac{a_{10} - b_{10}}{1 + e^{{{- k_{10}}{H(t)}} + c_{10}}} + b_{10}}},} & (7)\end{matrix}$

respectively. The parameters a₉, b₉, c₉, k₉, a₁₀, b₁₀, c₁₀, and k₁₀ arepositive constants with a₉>b₉ and a₁₀>b₁₀. The hepcidin level H(t) attime t is determined by equation (5).

3.2.2. Precursor Cells (Compartment 2).

Inflow into the precursor cell compartment comes from plasma. Oneoutflow of this compartment leaves to the erythrocytes compartment anddescribes the amount of hemoglobin carried from reticulocytes, whichenter the blood stream. Thus, the flow f₂₃(⋅) from compartment 2 tocompartment 3, is defined by the cells that leave the bone marrow andthe amount of hemoglobin these cells are carrying. These quantities canbe calculated from the population model and

k ₂₃ a ₂=ƒ₂₃(v ^(s)(E(t))s(t,μ _(max) ^(s),τ^(s)),h ^(s)(μ_(max) ^(s))),

where v^(s)(E(t))s(t, μ^(s) _(max), τ^(s)) determines how many cells areleaving the bone marrow at a certain time t and the distributionfunction h^(s)(μ^(s) _(max)) at the maximum maturity level μ^(s) _(max)indicates how much hemoglobin they carry. Here v^(s)(E(t)) is thematuration velocity of reticulocytes depending on EPO, s(t, μ^(s)_(max), τ^(s)) is the density of reticulocyte population with maximalmaturity μ^(s) _(max), and τ^(S)∈[τ^(s) _(min), τ^(s) _(max)] is thenumber of transferrin receptors at time t. The hemoglobin distributionfunction h^(s)(μ^(s) _(max)) arises as a result of the population model,in particular the precursor cell population part.

The second outflow of this compartment leaves to the macrophages and isdefined by the precursor cells that die. Hence, the flow fromcompartment 2 to compartment 4 is defined by the amount of iron thatdying cells carry, i.e.,

(8)

${{k_{24}a_{4}{\overset{\mu_{\max}^{r}}{\int\limits_{\mu_{\min}^{r}}}{\overset{h_{\max}^{r}}{\int\limits_{h_{\min}^{r}}}{\overset{\tau_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}}{h^{r}{\alpha^{r}\left( {\mu^{r},h^{r}} \right)}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}}} + {\overset{\mu_{\max}^{s}}{\int\limits_{\mu_{\min}^{s}}}{\overset{h_{\max}^{s}}{\int\limits_{h_{\min}^{s}}}{\underset{\tau_{\min}^{s}}{\int\limits^{\tau_{\max}^{s}}}{h^{s}{\alpha^{s}\left( {\mu^{s},h^{s}} \right)}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}},$

where r(t, μ^(r), τ^(r)) and s(t, μ^(s), τ^(s)) are the densities of thetwo precursor cell populations. Further, μ¹∈[μ^(i) _(min), μ^(i) _(max)]denotes the maturity, h^(i)∈[h^(i) _(min), h^(i) _(max)] denotes theamount of hemoglobin carried by a cell and τ^(i)∈[τ^(i) _(min), τ^(i)_(max)] denotes the number of transferrin receptors for I=r,s. Thefunction α^(i)(μ^(i), h^(i)), i=r,s describes the mortality rate ofprecursor cells, where α^(r)=α^(s) is determined by equations (3) and(4).

Altogether the change of amount of iron in the iron-precursorscompartment is given by

${\frac{d}{dt}{a_{2}(t)}} = {{{d_{0}\left( {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{\tau^{r}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{\tau^{s}{s\left( {t,{\mu^{s}h^{s}}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}} \right)}{a_{1}(t)}} - {f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s},\tau^{s}} \right)}},{h^{s}\left( \mu_{\max}^{s} \right)}} \right)} - {\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{a^{r}\left( {\mu^{r},h^{r}} \right)}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{\alpha^{s}\left( {\mu^{s},h^{s}} \right)}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d{\tau^{s}.}}}}}}$

Note that it is not necessary to determine the function f₂₃, becausethere is no need to explicitly observe the rate of change of theprecursor cell compartment. The amount of iron a₂(t) in this compartmentat time t can be directly computed from the population classes. It isdefined by the total number of precursor cells and the amount ofhemoglobin they carry:

${a_{2}(t)} = {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d{\tau^{s}.}}}}}}$

However, it is important to know the flux of hemoglobin to themacrophage compartment which arises from dying cells, because this flux(see equation (8)) is a source term in the macrophage compartment.

3.2.3. Erythrocytes (Compartment 3).

Inflow into the erythrocytes compartment comes from the precursor cellcompartment, and iron leaves the compartment to the macrophagescompartment (phagocytosis of senescent cells, random break-down of cellsand neocytolysis, shedding of hemoglobin containing vesicles)

${{k_{34}a_{3}} = {{\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}} + {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}}}},$

or is excreted via bleeding

${k_{3,{out}}a_{3}} = {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha_{bleed}^{m}(t)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{{dh}^{m}.}}}}$

Here α^(m) denotes the mortality rate of erythrocytes due tophagocytosis of senescent cells, random break-down of cells andneocytolysis. Furthermore, h^(m) is the amount of hemoglobin carried bya cell, m(t, μ^(m), h^(m)) is the population density of erythrocyteswith cell age μ^(m) and cellular hemoglobin h^(m) at time t and α^(m)_(bleed)(t) is the rate with which RBCs are lost due to bleeding. Thus,the change of amount of iron in the erythrocyte compartment is

${\frac{d}{dt}{a_{3}(t)}} = {{f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s},\tau^{s}} \right)}},{h^{s}\left( \mu_{\max}^{s} \right)}} \right)} - {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha_{bleed}^{m}(t)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}} - {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}} - {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{{dh}^{m}.}}}}}$

Similar to compartment 2, it is not necessary to observe the rate ofchange of the erythrocyte cell compartment to compute the amount of irona₃(t) at time t. The quantity a₃(t) can be directly computed from thepopulation class and is defined by the total number of circulating redblood cells and the amount of hemoglobin they carry:

${a_{3}(t)} = {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{m\left( {t,\mu^{m},h^{m}} \right)}{dh}^{m}d{\mu^{m}.}}}}$

However, the term k₃₄a₃(t) needs to be known explicitly (see equation(9)). This flux, which is the rate of hemoglobin released byerythrocytes, appears as a source term in the macrophage compartment.

3.2.4. Macrophages (Compartment 4).

Inflow into the macrophages compartment comes from the erythrocytescompartment (erythrocytes phagocytosed by macrophages, iron set free byrandom break-down of RBCs and shedding of hemoglobin containingvesicles) and the precursor cells compartment. Outflow of thiscompartment leaves to the plasma and to the storage compartment.Therefore, the change of iron in the compartment is

${{\frac{d}{dt}{a_{4}(t)}} = {{\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}} + {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}h^{m}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}} + {\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}},$

where k₄₁(H(t)) is given by equation (6).

3.2.5. Storage (Compartment 5.)

Inflow into the storage compartment comes from the plasma and themacrophages compartment and iron leaves this compartment to the plasmacompartment. Hence, the change of iron in the storage compartment is

${{\frac{d}{dt}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}},$

where k₅₁(H(t)) is given in equation (7).

The iron submodel can be applied independently of erythropoiesis to avariety of iron disorders, such as true iron deficiency states (e.g.,bleeding, malabsorption), functional iron deficiencies (e.g., anemias ofchronic disease, usually with chronic inflammatory disorders ormalignancies), or iron overload disorders (e.g., all forms ofhemochromatosis, iron intoxication, polytransfusions of blood, chronichepatopathies, and hyper plastic erythroid disorders such asthalassemias).

4. A Complete Listing of the Equations for the Mathematical Model.

The model consists of the following structured population equations:

${{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},{{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = \left( {{\beta^{q} - {{\alpha^{q}\left( {E(t)} \right)}{q\left( {t,\mu^{q}} \right)}}},} \right.}$${{\frac{\partial}{\partial t}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} + {{v_{2}^{r}\left( {\tau,{a_{1}(t)}} \right)}{\frac{\partial}{\partial h^{r}}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}}}},{+ {\frac{\partial}{\partial\tau^{r}}\left( {{{{v_{3}^{r}\left( {\mu^{r},h^{r},\tau^{r}} \right)}{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}} = {{8\beta{u\left( {t,\mu^{r},{2h^{r}},{2\tau^{r}}} \right)}} - {\left( {\beta + {\alpha^{r}\left( {\mu^{r},h^{r}} \right)}} \right){r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}}}},{{{\frac{\partial}{\partial t}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} + {{v_{1}^{s}\left( {E(t)} \right)}{\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{r},h^{r},\tau^{r}} \right)}}} + {{v_{2}^{s}\left( {\tau^{s},{a_{1}(t)}} \right)}{\frac{\partial}{\partial h^{s}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}} + {\frac{d}{d\tau^{s}}\left( {{v_{3}^{s}\left( {\mu^{s},h^{s},\tau^{s}} \right)}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}} \right)}} = {{- {\alpha^{s}\left( {\mu^{s},h^{s}} \right)}}{s\left( {t,\mu^{s},h^{s},\tau^{s}} \right)}}},{{{\frac{\partial}{\partial t}{m\left( {t,\mu^{m},h^{m}} \right)}} + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}} + {{v_{2}^{s}\left( \mu^{m} \right)}{\frac{\partial}{\partial h^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}}}} = {{- {\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},} \right.}}$${{\frac{d}{dt}{E^{end}(t)}} = {{\frac{1}{TBV}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},{{\frac{d}{dt}{E^{ex}(t)}} = {{\frac{1}{TBV}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{{E^{ex}(t)}.}}}}$

Further, the iron model can be described as follows:

${\frac{d}{dt}{a_{1}(t)}} = {k_{41}\left( {{{H(t)}a_{4}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} + {k_{51}\left( {{{{H(t)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{0}\left( {{\int\limits_{\mu_{\min}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{\tau^{r}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{\tau^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}} \right)}{{a_{1}(t)}.{a_{2}(t)}}}} = {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d{\tau^{s}.}}}}}}} \right.}} \right.}$${a_{3}(t)} = {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{m\left( {t,\mu^{m},h^{m}} \right)}{dh}^{m}d{\mu^{m}.}}}}$${{\frac{d}{dt}{a_{4}(t)}} = {{\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{h^{m}{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}} + {\int\limits_{\mu_{\min}^{m}}^{\mu_{\max}^{m}}{\int\limits_{h_{\min}^{m}}^{h_{\max}^{m}}{{v_{2}^{m}\left( \mu^{m} \right)}h^{m}{m\left( {t,\mu^{m},h^{m}} \right)}d\mu^{m}{dh}^{m}}}} + {\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{h^{r}{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{h^{s}{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}},$${{\frac{d}{dt}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}},{{\frac{d}{dt}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{{H(t)}.}}}}$

The boundary conditions are

${{{p\left( {t,0} \right)} = {S_{0}\left( {E(t)} \right)}},{{q\left( {t,\mu_{\min}^{q}} \right)} = {p\left( {t,\mu_{\min}^{p}} \right)}},{{r\left( {t,\mu_{\min}^{r},0,\tau^{r}} \right)} = {{q\left( {t,\mu_{\max}} \right)}{\eta\left( \tau^{r} \right)}}},{{r\left( {t,\mu^{r},h^{r},\tau^{r}} \right)} = 0},{{{if}\left( {\mu^{r},h^{r},\tau^{r}} \right){is}a{boundary}{point}{and}h^{r}} > 0},{{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\min}^{s},h^{s},\tau^{s}} \right)}} = {r\left( {t,\mu_{\max}^{r},h^{r},\tau^{r}} \right)}},{{m\left( {t,\mu^{m},h_{\min}^{m}} \right)} = 0}}{{m\left( {t,0,h^{m}} \right)} = \left\{ {\begin{matrix}{{v^{S}\left( {E(t)} \right)}{\int}_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{s\left( {t,\mu_{\max}^{s},h^{s},\tau^{s}} \right)}d\tau^{s}} & {{{{for}{E(t)}} \leq E^{*}},} \\{{\int}_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{r\left( {t,\mu_{\max}^{r},h^{r},\tau^{r}} \right)}d\tau^{r}} & {{{for}{E(t)}} > E^{*}}\end{matrix},} \right.}$

and the initial values are given by

p(0,μ^(p))=p ₀(μ^(p)),

q(0,μ^(q))=q ₀(μ^(q)),

r(0,μ^(r) ,h ^(r),τ^(r))=r ₀(μ^(r) ,h ^(r),τ^(r)),

s(0,μ^(s) ,h ^(s),τ^(s))=s ₀(μ^(s) ,h ^(s),τ^(s)),

m(0,μ^(m) ,h ^(m))=m ⁰(μ^(m) ,h ^(m)),

E ^(end)(0)=E ₀ ^(end),

E _(ex)(0)=E ₀ ^(ex),

a ₁(0)=a _(1,0),

a ₂(0)=a _(2,0),

a ₃(0)=a _(3,0),

a ₄(0)=a _(4,0),

a ₅(0)=a _(5,0),

H(0)=H ₀.

Moreover, there are number of auxiliary equations for the cellpopulation model

${{S_{0}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + e^{{{- k_{1}}{E(t)}} + c_{1}}} + b_{1}}},{{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + e^{{k_{2}{E(t)}} - c_{2}}} + b_{2}}},{{v_{2}^{r}\left( {\tau^{r},a_{1}} \right)} = {d_{2}\tau^{r}a_{1}}},{{v_{3}^{r}\left( {\mu^{r},h^{r},\tau^{r}} \right)} = {{a_{3}{\tau^{r}\left( {h_{*}^{r} - {b_{3}\left( h^{r} \right)}^{2}} \right)}} - \frac{c_{3}{\max\left( {\tau^{r},0} \right)}\left( {\mu^{r} - \mu_{\max}^{q}} \right)}{\left( {1 + h^{r}} \right)}}},{{\alpha^{r}\left( {\mu^{r},h^{r}} \right)} = \frac{{b_{4}\left( {h^{**} - h^{r}} \right)}^{2}}{1 + \left( {\mu^{**} - \left( {\mu^{r} - \mu_{\max}^{q}} \right)} \right)^{2}}},{{{\int}_{\tau_{\min}}^{\tau_{\max}}{\eta(\tau)}d\tau} = 1},{{v_{1}^{s}\left( {E(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + e^{{{- k_{5}}{E(t)}} + c_{5}}} + b_{5}}},{{v_{2}^{s}\left( {\tau^{s},a_{1}} \right)} = {d_{2}\tau^{s}a_{1}}},{{v_{3}^{s}\left( {\mu^{s},h^{s},\tau^{s}} \right)} = {{a_{3}{\tau^{s}\left( {h_{*}^{s} - {b_{3}\left( h^{s} \right)}^{2}} \right)}} - \frac{c_{3}{\max\left( {\tau^{s},0} \right)}\left( {\mu^{s} - \mu_{\max}^{r}} \right)}{\left( {1 + h^{s}} \right)}}},{{\alpha^{s}\left( {\mu^{s},h^{s}} \right)} = \frac{{b_{4}\left( {h^{**} - h^{s}} \right)}^{2}}{1 + \left( {\mu^{**} - \left( {\mu^{s} - \mu_{\max}^{r}} \right)} \right)^{2}}},{{v_{2}^{m}\left( \mu^{m} \right)} = {{b_{6}\left( {\mu^{m} - \mu_{*}^{m}} \right)}^{2} - c_{6}}},{{\alpha^{m}\left( {{E(t)},\mu^{m},h^{m}} \right)} = \left\{ {\begin{matrix}{\begin{matrix}{{\gamma^{m}\left( {\mu^{m},\mu^{m}} \right)} + {{\min\left( {\frac{c_{E}}{{E(t)}^{k_{E}}},b_{E}} \right)}{for}}} \\{{{E(t)} < \tau_{E}},{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}}}\end{matrix},} \\{{\gamma^{m}\left( {\mu^{m},\mu^{m}} \right)}{otherwise}}\end{matrix},{{\gamma^{m}\left( {\mu^{m},\mu^{m}} \right)} = \left\{ {\begin{matrix}{\begin{matrix}{{{\alpha_{r}^{m}{for}\mu^{m}} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},} \\{h^{m} \in \left\lbrack {{h_{\min}^{m} + \delta},h_{\max}^{m}} \right\rbrack}\end{matrix},} \\{\begin{matrix}{{{{\gamma_{\mu}^{m}\left( \mu^{m} \right)}{\gamma_{h}^{m}\left( h^{m} \right)}{for}\mu^{m}} \in \left( {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)},} \\{h^{m} \in \left( {h_{\min}^{m},{h_{\min}^{m} + \delta}} \right)}\end{matrix},} \\{{{\infty{for}\mu^{m}} \geq \mu_{\max}^{m}},{h^{m} \leq h_{\min}^{m}}}\end{matrix},{{\gamma_{\mu}^{m}\left( \mu^{m} \right)} = \frac{3\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}}},{{\gamma_{h}^{m}\left( h^{m} \right)} = \frac{3\delta^{2}}{\left( h^{m} \right)^{2} - {2\left( {h_{\min}^{m} - \delta} \right)h^{m}} + {\left( {h_{\min}^{m} - {2\delta}} \right)h_{\max}^{m}}}},{{E_{in}^{end}(t)} = {\frac{a_{7} - b_{7}}{a + e^{\frac{k_{7}{{Hb}(t)}}{({TBV})} - c_{7}}} + b_{7}}},{{{Hb}(t)} = {{\int}_{h_{\min}^{m}}^{h_{\max}^{m}}{m\left( {t,\mu^{m},h^{m}} \right)}{dh}^{m}}},{{E(t)} = {{E^{ex}(t)} + {E^{in}(t)}}},} \right.}} \right.}$

and also some auxiliary equations for the iron compartment model

${{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{8} - b_{8}}{1 + e^{{{- k_{8}}{H(t)}} + c_{8}}} + b_{8}}},{{H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + e^{{{- k_{H}}{C(t)}} + c_{H}}} + b_{H}}},{{C(t)} = {{{- c_{U}}{U(t)}} + {c_{a_{5}}{a_{5}(t)}} + {c_{i}{i(t)}}}},$${{U(t)} = {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{\int\limits_{h_{\min}^{r}}^{h_{\max}^{r}}{\int\limits_{\tau_{\min}^{r}}^{\tau_{\max}^{r}}{{r\left( {t,\mu^{r},h^{r}} \right)}d\mu^{r}{dh}^{r}d\tau^{r}}}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{\int\limits_{h_{\min}^{s}}^{h_{\max}^{s}}{\int\limits_{\tau_{\min}^{s}}^{\tau_{\max}^{s}}{{s\left( {t,\mu^{s},h^{s}} \right)}d\mu^{s}{dh}^{s}d\tau^{s}}}}}}},$${{k_{41}\left( {H(t)} \right)} = {\frac{a_{9} - b_{9}}{1 + e^{{{- k_{9}}{H(t)}} + c_{9}}} + b_{9}}},{{{and}{k_{51}\left( {H(t)} \right)}} = {\frac{a_{10} - b_{10}}{1 + e^{{{- k_{10}}{H(t)}} + c_{10}}} + {b_{10}.}}}$

Exemplification of the Above Described Models

5. Simulations for the Iron Model.

TABLE 3 Iron compartments in normal man (estimates for an average man:Height: 177 cm, weight: 70 kg. See Williams Hematology Table 40.1)Compartment Iron content (mg) Total body iron (%) Hemoglobin iron 220067 (Precursors + Erythrocytes) Storage iron 1000 27 Labile pool(macrophages) 80 2.2 Transport iron (plasma) 3 0.08 Other tissue iron138 3.72

TABLE 4 Parameter values for the iron model in a steady state. ParameterMeaning Value k₁₂ transfer rate from plasma to precursor cells 8.5 k₁₅transfer rate from plasma to storage 1.62 k₂₃ transfer rate fromprecursor cells to erythrocytes 0.1199 k₂₄ transfer rate from precursorcells to macrophages 0.0101 k₃₄ transfer rate from erythrocytes tomacrophages 0.0012 k₄₁ transfer rate from macrophages to plasma 0.29965k₄₅ transfer rate from macrophages to plasma 0.0256 k₅₁ transfer ratefrom storage to plasma 0.007

Values for the parameters of the iron submodel are assigned as describedbelow. First, parameters are formed that are able to reflect thesituation in a steady state. For this purpose, the population model andthe iron model are decoded and the compartment model describing ironhomeostasis is the only model investigated. Further, because of steadystate, all transfer rates are assumed to be constant. Then, values areassigned to the transfer rates, such that the flow of iron from onecompartment to another, and also the total amount of iron in eachcompartment, are assumed to be within physiological ranges. Moreover,for instance in Williams Hematology, a table can be found where the ironcontent for a 70 kg man is presented. See Table 3 for a list of thevalues.

Parameter values chosen based on this information are shown in Table 4.Note that in the steady state, losses of iron and absorption of iron areequal and, therefore, are not considered here. Thus, for the moment,iron homeostasis is considered as a closed system. Using these parametervalues, a simulation was started close to the steady state values foreach compartment. In addition, a situation was considered where moreiron is incorporated in hemoglobin and, as a consequence, storage ironis decreased. Since, in this case the model is run with constanttransfer rates and no feedback or regulation is taken into account, ittakes quite a long time until the system finally reaches its steadystate.

Further Development of the Erythropoiesis and Iron Hemeostasis Models

The models presented in Sections 2 and 3 above are very close tophysiology and incorporate a lot of physiological mechanisms involved inerythropoiesis and iron homeostasis. As a consequence, the models arevery extensive and lead to a class of extremely complicated mathematicalequations. To solve the models, special numerical schemes have beendeveloped to find numerical approximations for the population equations.See D. H. Fuertinger. A model for erythropoiesis. PhD thesis, Universityof Graz, Austria, 2012 (hereinafter “Fuertinger Thesis”). Although theapproach chosen is very efficient, it is extremely computationallyintensive, because systems of several thousands ordinary differentialequations (ODE) have to be solved at every time step. Therefore, at thepresent time, simulations are only realizable on high performancecomputers, and thus not really practicable or viable in every dayclinical practice and patient care, where such computers are nottypically available. Since the iron component is such an important partof treatment of the anemia of renal failure, a simplified version of themodels presented in Sections 2 and 3 is described hereinafter. The modelpresented hereinafter is still comprehensive and demanding with regardto computational costs, although the complexity of the numericalapproximations reduces significantly. Instead of solving systems ofseveral thousand ODEs, less than about a hundred ODEs have to be solvedat every time step, which reduces the computational effort tremendously.

The first step is to restrict the population model presented in Section2 to population equations with only one structural attribute, namelycell age (also called maturity of the cell). Avoidance of the structuralattributes hemoglobin and transferrin receptors makes it necessary totake a step backward and reconsider how an iron model can be linked to amodel describing the erythroid cell populations. Several of theassumptions made in Sections 2 and 3 have to be revisited. It is stillimportant to keep track of how much iron is taken up by precursor cellsand how much hemoglobin is contained in every cell. The lack of thestructural feature hemoglobin, which allowed one to compute thisquantity relatively easily, requires that other mechanisms have to bedeveloped to replace it. The equations and assumptions are chosen suchthat the simplified model is as close to the one previously shown aspossible. Nevertheless, these ‘workarounds’ only enable the descriptionof the phenomenon of the accumulation of hemoglobin in cells and not somuch the physiological mechanisms involved. Whereas this might not makea difference in most cases, the more extensive model might be needed toexplain the development of RBC comprehensively for some patients withsevere changes in these mechanisms.

To ensure that the description below can be read independently, a fewimportant aspects of erythropoiesis and iron homeostasis are repeatedhereinafter.

6. Erythropoiesis Model.

The aim is to develop a mathematical model which incorporates the basicmechanisms governing erythropoiesis including the control actions oferythropoietin and the potentially restricting effect of iron. Thus, themodel has to include the regulatory processes for iron homeostasis, inparticular describing those which are connected to erythropoiesis. Theresult is a comprehensive mathematical model which is able to predictred blood cell mass under EPO- and iron-therapy.

As a first step, in this chapter a model is developed which focusessolely on the effects of erythropoietin on erythroid cells. Hence,throughout this section, erythropoiesis is assumed to be sufficientlysupplied with iron. Hereinafter, this model will be referred to as theErythropoiesis Model, described in PCT Application No.PCT/US/2012/054264 published as WO 2013/036836 A2 on Mar. 14, 2013,whose contents and teachings are incorporated herein by reference intheir entirety. The mathematical model developed in this section isapplicable for a number of different situations, as long as it isreasonable to assume that erythropoiesis is not impaired because of alack of iron availability. For instance, the Erythropoiesis Model isable to describe the recovery of red cell mass after blood donation, thereaction of the body to presurgical administration of ESAs, and changesin the number of erythrocytes of high altitude dwellers descending tosea level. See Fuertinger Thesis; D. H. Fuertinger, F. Kappel, S.Thijssen, N. W. Levin, and P. Kotanko. Journal of Mathematical Biology,66(6):1209-1240, 2013. However, the Erythropoiesis Model is onlyapplicable for a subpopulation of dialysis patients. 80-90% of dialysispatients treated with ESA suffer from functional iron deficiency at somestage in their therapy. See R. M. Schaefer and L. Schaefer. NephrologyDialysis Transplantation, 13:9-12, 1998. It has to be consideredcarefully whether, even when treated with iron compounds, it isreasonable to assume that iron supply is sufficient in patientssuffering from chronic kidney disease (CKD patients on dialysistreatment or not). Moreover, the Erythropoiesis Model helps tounderstand what the most important dynamics considering red blood cellproduction are. In Section 7 an Iron Model (compartment model) isdeveloped and in Section 7.2 it is explained how this model is linked tothe Erythropoiesis Model and extends it to a even more extensive anddetailed model for erythropoiesis, thus providing a solid basis for amore general description of the altered processes during the verycomplex pathological situation of anemia of renal failure.

The Erythropoiesis Model is based on structured population modelsdescribing the different erythroid cell stages. Five differentpopulation classes of cells are considered. BFU-Es, CFU-Es,erythroblasts, marrow reticulocytes and mature erythrocytes circulatingin the bloodstream (including peripheral reticulocytes). In theErythropoiesis Model individual cells are distinguished according totheir maturity, which can also be referred to as cell age. Thecommitment to the erythroid lineage is an irreversible event. Adifferentiated cell cannot regress or switch into anotherdifferentiation pathway. Thus, once a multipotent stem cell is committedto the erythroid lineage, it undergoes the complete series ofdifferentiations until it becomes a mature red blood cell, or it dieseventually during this process. While maturing, the cell divides anumber of times. Hence, age-structured population models of the form

$\begin{matrix}{{\left. {{{\frac{\partial}{\partial t}{u\left( {t,\mu} \right)}} + {{v\left( {E(t)} \right)}{\frac{\partial}{\partial\mu}{u\left( {t,\mu} \right)}}}} = {\left( {{\beta\left( {E(t)} \right)},\mu} \right) - {\alpha\left( {{E(t)},\mu} \right)}}} \right){u\left( {t,\mu} \right)}},} & (10)\end{matrix}$

are used in order to describe the development of the cell populations.Here u(t, μ) denotes the population density of the cell population attime t with maturity μ. Further, the functions β and α describe thedivision rate and the rate of apoptosis of the cells, respectively. Thefunctions α and β a priori depend on the maturity μ and theconcentration of EPO E(t), respectively, at time t. The function vdescribes the maturation velocity and depends on the concentration E(t).

For the different population classes, the characteristic properties(proliferation rate, rate of apoptosis and maturation velocity) differ.While the cell matures, it changes its morphological characteristics,such as the number of EPO- and transferrin-receptors expressed on thesurface.

In the following sections, the assumptions that were made for thedifferent population classes are listed, briefly described, and themathematical equations that arise in consequence are stated. FIG. 9Agives an overview of the design and organization of the model.

6.1. Progenitor Cells: BFU-E and CFU-E Cells

Assumptions:

-   -   1. The number of cells, which commit to the erythroid lineage,        is constant.    -   2. Cells normally stay in this stage for 8 days (3 days BFU-E        and 5 days CFU-E).    -   3. EPO has no effect on the number of divisions or the rate of        apoptosis of BFU-E.    -   4. The proliferation rate of CFU-E cells is constant.    -   5. The rate of apoptosis of CFU-E depends highly on EPO levels.    -   6. The maturation velocities of BFU-E and CFU-E cells are        constant.

The process by which stem cells are recruited into proliferatingprogenitor population remains unclear. There are several hypothesesincluding an environmental dependency, that it is a random event, etc.For the moment, the number of stem cells entering the erythroid lineageis assumed to be independent of EPO and thus constant. The change inpopulation of the progenitor cells over time are described consideringtwo different classes of cells: namely BFU-E and CFU-E cells.

The earliest identifiable erythroid progenitor cell is the Burst-formingUnit Erythroid (BFU-E). At first these cells express only a very smallnumber of EPO receptors on the surface. See Wintrobe's ClinicalHematology. Thus, it is reasonable to assume EPO has no effect onproliferation or apoptosis of these cells. See Wu et al. In culture itlasts around 6-7 days until human BFU-E have all the functionalcharacteristics of the next cell stage—Colony-forming Unit Erythroidcells (CFU-E cells) (Assumption 2). See Williams Hematology.Morphologically, it is difficult to distinguish between those two typesof cells, because there are cells in between these two developmentalstages which show characteristic properties between BFU-E and CFU-Ecells. Therefore, a distinction is valid but artificial.

The BFU-E cell class is described by the following population equation

$\left\{ {\begin{matrix}{{{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},} \\{{{p\left( {t,0} \right)} = S_{0}},} \\{{p\left( {0,\mu^{p}} \right)} = {p_{0}\left( \mu^{p} \right)}}\end{matrix},} \right.$

where p(t, μ^(p)) is the population density of the cell class at time twith maturity μ^(p), 0≤μ^(p)≤μ^(p) _(max)=7, t>0. Further, β^(p) is aconstant proliferation rate and α^(p)≡0 (Assumption 3), S₀ describes thenumber of cells committing to the erythroid lineage (Assumption 1) andp₀(μ^(p)) is the population density at t=0.

Once a cell reaches the maximum age for BFU-E cells, it leaves thispopulation class and enters the CFU-E class. Consequently, there is acontinual flux of cells from one population class to the next one. CFU-Eare more rapidly dividing cells than BFU-E. See Wu et al. During thisstage, cells are very sensitive to EPO levels, and, under normalconditions, large numbers of generated CFU-E are not surviving. SeeWintrobe's Clinical Hematology. CFU-E are highly dependent on EPO toprevent them from apoptosis, i.e., the mortality for this populationclass depends on the EPO concentration. Altogether, the followingequations for the second class are obtained:

$\left\{ {\begin{matrix}\begin{matrix}{{{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},} \\{{{q\left( {t,\mu_{\min}^{q}} \right)} = {p\left( {t,\mu_{\max}^{p}} \right)}},}\end{matrix} \\{{q\left( {0,\mu^{q}} \right)} = {q_{0}\left( \mu^{q} \right)}}\end{matrix},} \right.$

where q(t, μ^(q)) is the population density of the CFU-E class at time twith maturity μ^(q), t>0 and 7=μ^(q) _(min)≤μ^(q)≤μ^(q) _(max)=13.Further on, β^(q) stands for a constant proliferation rate (Assumption4), α^(q) (E(t)) denotes the apoptosis rate depending on theEPO-concentration (Assumption 5), q(t, μ^(q) _(min))=p(t, μ^(p) _(max))describes the number of cells leaving the BFU-E cell stage and enteringthe CFU-E cell stage and q₀(μ^(q)) is the population density at t=0.

A sigmoid function is used to describe the rate of apoptosis,

${{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + e^{{k_{1}{E(t)}} - c_{1}}} + b_{1}}},$

where E(t) is the EPO concentration at time t and a₁, b₁, c₁ and k₁ arepositive constants with a₁>b₁. The function α^(q) monotonicallydecreases with increasing EPO concentration. Thus, a higher level of EPOcauses more cells to survive.

Note that, because of Assumption 6, the cell age of progenitor cellscoincides with the actual age of the cell, i.e., one definesv^(p)=v^(q)≡1, ∀t≥0.

6.2. Precursor Cells: Erythroblasts and Marrow Reticulocytes

Assumptions:

-   -   7. Cells stay in this stage for 6-8 days (5 days erythroblasts        and 1-3 days marrow reticulocytes).    -   8. The class erythroblasts consists of all cell stages from        proerythroblast to orthochromatophilic erythroblast.    -   9. EPO has no effect on the number of divisions or the rate of        apoptosis of erythroblasts. The proliferation rate of        erythroblasts is assumed to be constant.    -   10. The maturation velocity of erythroblasts is constant.    -   11. Reticulocytes mature but do not proliferate.    -   12. The maturation velocity of reticulocytes depends on EPO.    -   13. A constant portion of marrow reticulocytes is phagocytosed.

After a CFU-E differentiates to a proerythroblast, it takes aboutanother 6-8 days until the cell is released from the bone marrow intothe bloodstream. See Jandl. The various stages of maturation fromproerythroblast to orthochromatophilic erythroblast are referred to aserythroblasts. The cells undergo several mitotic divisions until, at thestage of orthochromatophilic erythroblast, they lose their ability todivide and enter a maturation period. The erythroblastic pyramids appearnormal, with no evidence of additional mitotic divisions, whenproduction increases, i.e., the proliferation of erythroblasts isassumed to be independent of EPO levels and defined to be constant. SeeWilliams Hematology.

Hence, the erythroblast class can be described by the following equation

$\left\{ {\begin{matrix}\begin{matrix}{{{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r}} \right)}}} = {\beta^{r}{r\left( {t,\mu^{r}} \right)}}},} \\{{{r\left( {t,\mu_{\min}^{r}} \right)} = {q\left( {t,\mu_{\min}^{q}} \right)}},}\end{matrix} \\{{r\left( {0,\mu^{r}} \right)} = {r_{0}\left( \mu^{r} \right)}}\end{matrix},} \right.$

where r(t, μ^(r)) is the population density of the erythroblasts at timet with maturity μ^(r), t>0 and 13=μ^(r) _(min)≤μ^(r)≤μ^(r) _(max)=18.Further, β^(r) is a constant proliferation rate and α^(s)≡0 (Assumption9), the maturation velocity v^(r)≡1 (Assumption 10), r(t, μ^(r)_(min))=q(t, μ^(q) _(max)) describes the number of cells leaving theCFU-E stage and entering the erythroblasts cell stage, and r₀(μ^(r)) isthe population density at t=0.

The differentiating process from orthochromatophilic erythroblasts tomarrow reticulocytes involves the extrusion of the cell nucleus.Reticulocytes are not capable of cell divisions (Assumption 11), i.e.,β^(s)≡0. In the Erythropoiesis Model, one does not account for animpaired erythropoiesis due to iron deficiency. Still, even when theerythroid cells in the bone marrow are sufficiently supplied with iron,not all precursor cells survive. Some of the reticulocytes die beforethey are released to the blood stream. See G. Barosi, M. Cazzola, C.Berzuini, S. Quaglini, and M. Stefanelli. British Journal ofHaematology, 61:357-370, 1985; M. Stefanelli, D. P. Bentley, I. Cavill,and H. P. Roeser. The American Journal of Physiology, 247:842-849, 1984.This is why Assumption 13 is made.

A raised EPO concentration shortens the marrow transit time of precursorcells. If the EPO level is elevated, marrow reticulocytes are releasedprematurely. This aspect is accounted for by allowing the maturationvelocity of reticulocytes to vary depending on erythropoietinconcentration. Thus, although the maximum cell age of marrowreticulocytes is fixed, the actual transit time for the cells variesbetween 1-3 days. Hence, altogether the transit time for precursor cellsis between 6-8 days (Assumption 7).

The population equation referring to marrow reticulocytes is

$\left\{ {\begin{matrix}\begin{matrix}{{{{\frac{\partial}{\partial t}{s\left( {t,\mu^{s}} \right)}} + {{v^{s}\left( {E(t)} \right)}{\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{s}} \right)}}}} = {{- \alpha^{s}}{s\left( {t,\mu^{s}} \right)}}},} \\{{{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\min}^{s}} \right)}} = {r\left( {t,\mu_{\max}^{r}} \right)}},}\end{matrix} \\{{s\left( {0,\mu^{s}} \right)} = {r_{0}\left( \mu^{s} \right)}}\end{matrix},} \right.$

where s(t, μ^(s)) is the reticulocytes population density at time t withmaturity μ^(s), t>0 and 18=μ^(s) _(min)≤μ^(s)≤μ^(s) _(max)=20. Furtheron, v^(s)(E(t)) is the maturation velocity depending on EPO (Assumption12), α^(s) denotes the rate with which reticulocytes are phagocytosed(Assumption 13), v^(s)(E(t))s(t, μ^(s) _(min))=r(t, μ^(r) _(max))describes the number of cells leaving the erythroblast cell stage andentering the reticulocyte cell stage, and r₀(μ^(s)) is the populationdensity at t=0.

A sigmoid function is used to describe the changes in maturationvelocity,

${{v^{s}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + e^{{{- k_{2}}{E(t)}} + c_{2}}} + b_{2}}},$

where E(t) is the EPO concentration at time t and a₂, b₂, c₂ and k₂ arepositive constants with a₂>b₂. The maturation velocity increases with arising concentration of EPO. Note that a slower maturation velocitycauses the cells to reach the maximum cell age μ^(s) _(max) at a laterpoint, whereas a faster maturation velocity shortens the transit time,i.e., the cells reach μ^(s) _(max) earlier.

6.3. Erythrocytes.

Assumptions:

-   -   14. Erythrocytes and blood reticulocytes are subsumed in one        class.    -   15. Cells mature but do not proliferate.    -   16. There is a fixed random daily break-down of red blood cells        (not to be confused with loss of erythrocytes due to        senescence).    -   17. A drop in EPO concentration beneath a threshold precipitates        neocytolysis.    -   18. Erythrocytes with age between 14-21 days are likely to be        affected by neocytolysis.    -   19. Cells are phagocytosed when they reach the maximum age.    -   20. The maximum age of erythrocytes in a healthy persons is 120        days.

The reticulocytes are released from the bone marrow into the bloodstream and within 1-2 days they mature to erythrocytes. Circulating redblood cells have no nuclei and therefore they are not able toproliferate and they are not able to repair themselves. Thus, the lifespan of these cells is limited. In healthy adults the average life spanis about 120 days, but it can significantly shorten in some pathologies.See e.g., Jandl. If not otherwise stated, 120 days is used as themaximal age of erythrocytes for the Erythropoiesis Model.

Cells can be lost for different reasons: due to internal or externalbleeding, because of random daily-breakdown, because neocytolysis istriggered, and, because of eryptosis of senescent cells.

Altogether, the following equation is obtained for red blood cells:

$\left\{ {\begin{matrix}\begin{matrix}{{{{\frac{\partial}{\partial t}{m\left( {t,\mu^{m}} \right)}} + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m}} \right)}}} = {{- {\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},} \\{{{m\left( {t,0} \right)} = {{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s}} \right)}}},}\end{matrix} \\{{m\left( {0,\mu^{m}} \right)} = {m_{0}\left( \mu^{m} \right)}}\end{matrix},} \right.$

where m(t, μ^(m)) is the population density for the erythrocyte class attime t with maturity μ^(m), t>0 and 0=μ^(m) _(min)≤μ^(m)≤μ^(m)_(max)=120 (Assumption 20). Moreover, m(t,0)=v^(s)(E(t))s(t, μ^(s)_(max)) describes the number of reticulocytes entering the blood streamand m₀(μ^(m)) is the population density at t=0. In the populationequation α^(m)(E(t), μ^(m)) denotes a random daily break-down andneocytolysis

${\alpha^{m}\left( {{E(t)},\mu^{m}} \right)} = \left\{ {\begin{matrix}{{{{\gamma^{m}\left( \mu^{m} \right)} + {{\min\left( {\frac{c_{E}}{{E(t)}^{k_{E}}},b_{E}} \right)}{for}{E(t)}}} < \tau_{E}},{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}},} \\{{\gamma^{m}\left( \mu^{m} \right)}{otherwise}}\end{matrix},} \right.$

where b_(E), c_(E), k_(E) are positive constants, τ_(E) is the thresholdbeneath which neocytolysis is triggered (Assumption 17, and [μ_(min)^(m,n),μ_(max) ^(m,n)]=[14,21](Assumption 18) the age-interval duringwhich cells are affected by neocytolysis, and t≥0. Additionally, it ispossible to add a further term α^(m) _(bleed)(t) to the mortality ratein case of bleeding. Since it is assumed that a cell is phagocytosedwhen it reaches its maximum age (Assumption 19), the mortality rateγ^(m)(μ^(m)) needs to be chosen such that γ^(m)(μ^(m))=α^(m) _(rand) forμ^(m) _(min)≤μ^(m)≤μ^(m) _(max)−δ, with δ>0 sufficiently small and ∫_(μ)_(max) _(m) _(−δ) ^(μ) ^(max) ^(m) γ^(m)(μ^(m))=∞. Here α^(r) _(rand) isa random daily break-down (Assumption 16). A possible choice forγ^(m)(μ^(m)) is

${\gamma^{m}\left( \mu^{m} \right)} = \left\{ {\begin{matrix}{{{\alpha_{rand}^{m}{for}\mu^{m}} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},} \\{{{\frac{3\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}}{for}\mu^{m}} \in \left( {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)},} \\{{\infty{for}\mu^{m}} \geq \mu_{\max}^{m}}\end{matrix}.} \right.$

6.4. Erythropoietin

Assumptions:

-   -   21. Release of EPO is controlled by a negative feedback        mechanism according to the oxygen content.    -   22. Oxygen carrying capacity is directly proportional to the        number of erythrocytes.    -   23. The degradation rate of EPO is constant.    -   24. There is a slight delay in reaction of the EPO production        rate to the number of RBC but this is negligible compared to the        duration of development of erythrocytes.

The kidneys adjust the release of EPO according to the oxygen carryingcapacity of the blood. If blood oxygen content is lower than normal,then EPO production increases and vice versa. Thus, the production ofEPO is controlled by a negative feedback mechanism and allows for morered blood cells to be developed in case of an undersupply of the bodywith oxygen by opposing programmed cell death of erythroid progenitorcells. Further, if EPO decreases beneath a certain threshold,neocytolysis is triggered, that is, a process wherein macrophages startto phagocytose young erythrocytes (neocytes).

A sigmoid function which depends on blood oxygen partial pressure isused to model the feedback involving the release of erythropoietinE_(in) ^(end)(t) from the kidneys into plasma. As a consequence ofAssumption 22, the amount of E_(in) ^(end)(t) of EPO released by thekidney per unit time can be directly computed by use of the totalpopulation of erythrocytes M(t)=∫₀ ^(μ) ^(max) ^(m) (t, μ^(m))dμ^(m).Recall that this class consists of all circulating red blood cells(Assumption 14).

${{E_{in}^{end}(t)} = {\frac{a_{3} - b_{3}}{1 + e^{{k_{3}{\overset{\sim}{M}(t)}} - c_{3}}} + b_{3}}},$

where

${\overset{\sim}{M}(t)} = \frac{10^{- 8}{M(t)}}{TBV}$

is a erythrocytes “concentration”. Here, TBV is the total blood volume.The constants a₃, b₃, c₃, k₃ are positive and satisfy a₃>b₃. Thefunction E_(in) ^(end)(t) is monotonically decreasing. Thus, the releaseof EPO increases if the number of circulating red blood cells decreases(Assumption 21). The dynamics of the endogenous EPO concentrationE^(end)(t) in plasma are described by the following ordinarydifferential equation:

$\begin{matrix}{{{\frac{d}{dt}{E^{end}(t)}} = {{\frac{1}{TBV}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},} & (11)\end{matrix}$

where E^(end)(t) is the endogenous EPO concentration in plasma, E_(in)^(end) is the amount of EPO released by the kidneys and c_(deg) ^(ex)describes the constant degradation rate of endogenous EPO (Assumption23).

The degradation rate c_(deg) ^(ex) for exogenous EPO differs from theone for endogenous EPO and varies according to the kind of ESAadministered. Therefore, an additional ODE is needed to describe thechange in the plasma concentration of an ESA

$\begin{matrix}{{\frac{d}{dt}{E^{ex}(t)}} = {{\frac{1}{TBV}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{{E^{ex}(t)}.}}}} & (12)\end{matrix}$

Here E_(in) ^(ex)(t) is the rate at which the artificial hormone isadministered and c_(deg) ^(ex) is the rate with which the exogenoushormone is degraded. In intravenous administration, the total amount ofthe agent is injected into a vein, within a very short time interval. Inthis case E_(in) ^(ex)(t) can be approximated by E₀ ^(ex)(t)δ_(t0)(t),where E₀ ^(ex) is the amount of artificial hormone administered andδ_(t0)(t) is the Dirac delta impulse located at to, the time when theadministration takes place. The overall concentration of EPO in bloodconsists of the naturally produced erythropoietin in the body and theadministered ESA

E(t)=E ^(ex)(t)+E ^(end)(t)

6.5. Mathematical Model Erythropoiesis Model

For the convenience of the reader, all equations of the ErythropoiesisModel are collected in this section:

${{{\frac{\partial}{\partial t}{p\left( {t,\mu^{p}} \right)}} + {\frac{\partial}{\partial\mu^{p}}{p\left( {t,\mu^{p}} \right)}}} = {\beta^{p}{p\left( {t,\mu^{p}} \right)}}},$${{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},$${{{\frac{\partial}{\partial t}{r\left( {t,\mu^{r}} \right)}} + {\frac{\partial}{\partial\mu^{r}}{r\left( {t,\mu^{r}} \right)}}} = {\beta^{r}{r\left( {t,\mu^{r}} \right)}}},$${{{\frac{\partial}{\partial t}{s\left( {t,\mu^{s}} \right)}} + {{v^{s}\left( {E(t)} \right)}{\frac{\partial}{\partial\mu^{s}}{s\left( {t,\mu^{s}} \right)}}}} = {{- \alpha^{s}}{s\left( {t,\mu^{s}} \right)}}},$${{{\frac{\partial}{\partial t}{q\left( {t,\mu^{q}} \right)}} + {\frac{\partial}{\partial\mu^{q}}{q\left( {t,\mu^{q}} \right)}}} = {\left( {\beta^{q} - {\alpha^{q}\left( {E(t)} \right)}} \right){q\left( {t,\mu^{q}} \right)}}},$${{{\frac{\partial}{\partial t}{m\left( {t,\mu^{m}} \right)}} + {\frac{\partial}{\partial\mu^{m}}{m\left( {t,\mu^{m}} \right)}}} = {{- {\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}}{m\left( {t,\mu^{m}} \right)}}},$${{\frac{d}{dt}{E^{end}(t)}} = {{\frac{1}{TBV}{E_{in}^{end}(t)}} - {c_{\deg}^{end}{E^{end}(t)}}}},$${\frac{d}{dt}{E^{ex}(t)}} = {{\frac{1}{TBV}{E_{in}^{ex}(t)}} - {c_{\deg}^{ex}{E^{ex}(t)}{where}}}$${{\alpha^{q}\left( {E(t)} \right)} = {\frac{a_{1} - b_{1}}{1 + e^{{k_{1}{E(t)}} - c_{1}}} + b_{1}}},$${v^{s}\left( {E(t)} \right)} = {\frac{a_{2} - b_{2}}{1 + e^{{{- k_{2}}{E(t)}} + c_{2}}} + {b_{2}.{and}}}$${\alpha^{m}\left( {{E(t)},\mu^{m}} \right)} = \left\{ {\begin{matrix}{{{{\gamma^{m}\left( \mu^{m} \right)} + {{\min\left( {\frac{c_{E}}{{E(t)}^{k_{E}}},b_{E}} \right)}{for}{E(t)}}} < \tau_{E}},{\mu_{\min}^{m,n} \leq \mu^{m} \leq \mu_{\max}^{m,n}},} \\{{\gamma^{m}\left( \mu^{m} \right)}{otherwise}}\end{matrix},{with}} \right.$${\gamma^{m}\left( \mu^{m} \right)} = \left\{ {\begin{matrix}{{{\alpha_{rand}^{m}{for}\mu^{m}} \in \left\lbrack {\mu_{\min}^{m},{\mu_{\max}^{m} - \delta}} \right\rbrack},} \\{{{\frac{3\alpha_{rand}^{m}\delta^{2}}{\left( \mu^{m} \right)^{2} - {2\left( {\mu_{\max}^{m} + \delta} \right)\mu^{m}} + {\left( {\mu_{\max}^{m} + {2\delta}} \right)\mu_{\max}^{m}}}{for}\mu^{m}} \in \left( {{\mu_{\max}^{m} - \delta},\mu_{\max}^{m}} \right)},} \\{{\infty{for}\mu^{m}} \geq \mu_{\max}^{m}}\end{matrix}.} \right.$

Further, the boundary conditions are given by

p(t,0)=S ₀,

q(t,μ _(min) ^(q))=p(t,μ _(max) ^(p)),

r(t,μ _(min) ^(r))=q(t,μ _(max) ^(q)),

v ^(s)(E(t))s(t,μ _(min) ^(s))=r(t,μ _(max) ^(r)),

m(t,0)=v ^(s)(E(t))s(t,μ _(max) ^(s)),

m(0,μ^(m))=m ₀(μ^(m)),

and the initial values are

p(0,μ^(p))=p ₀(μ^(p)),

q(0,μ^(q))=q ₀(μ^(q)),

r(0,μ^(r))=r ₀(μ^(r)),

s(0,μ^(s))=r ₀(μ^(s)),

E ^(end)(0)=E ₀ ^(end),

E _(ex)(0)=E ₀ ^(ex).

The feedback is described by

${{E_{in}^{end}(t)} = {\frac{a_{3} - b_{3}}{1 + e^{{k_{3}{\overset{\sim}{M}(t)}} - c_{3}}} + b_{3}}},$

and one has

E(t)=E ^(ext)(t)+E ^(end)(t).

7. Iron Homeostasis.

As described above, much of the iron in the human body can be found incirculating erythrocytes incorporated in hemoglobin. However, most ofthe iron used for the production of new erythrocytes comes fromhemoglobin recycling and not from absorption. When a senescent red cellis phagocytosed, the macrophage internalizes hemoglobin, degrades it andexports it back out into blood circulation.

Iron is stored in form of ferritin or hemosiderin, where the latter ismore a kind of a long-term storage, because iron within deposits ofhemosiderin is very poorly available. Storage sites can be found mainlyin the liver, the spleen and the bone marrow but, for instance, smallamounts of hemosiderin can be found in almost any tissue. There aredifferences between men and women in the amount of iron that is stored(Females about 200-400 mg, Males about 1000 mg). See Crepaldi.

Iron is carried around in plasma by a transport protein, which is calledtransferrin. Transferrin is produced by the liver. Cells that are inneed for iron express transferrin receptors, and after transferrin hasformed a complex with this receptor, iron is transported into the cell.Transferrin exists in three different forms: iron-free (apotransferrin),one iron (monoferric transferrin) and two iron (differic transferrin)atoms bound. The concentration of iron and the amount of transferrinpresent in plasma determine the relative abundance of each form.Although iron bound to transferrin amounts to only about 3 mg of totalbody iron, it is the most important iron pool with the highest rate ofturnover, about 30 mg/day. Approximately 80% of this iron is transportedto the bone marrow and used by erythroid cells. See Lieu et al. Theamount of transferrin-bound iron is determined by three processes:macrophage iron recycling, (hepatic) iron storage and duodenal ironabsorption. See Wrighting.

In the last 15 years, understanding of the regulation of systemic ironhomeostasis has changed substantially. The antimicrobial peptidehepcidin, which is secreted by the liver, was identified to be thesystemic iron regulator. Hepcidin inhibits iron transport by binding tothe iron exporter ferroportin. Hence, it keeps gut enterocytes fromsecreting iron into circulation, thereby functionally reducing ironabsorption. Further, it prevents iron release from macrophages andblocks the stores as well, by shutting off the means of transport out ofcells.

There exist some sensing mechanisms in the liver that gauge the amountof iron needed to support erythropoiesis. Note that hepatocytes in theliver, which secrete hepcidin, also release transferrin into circulationand are one of the major iron storage sites. If the concentration ofiron stores are too low and/or erythropoiesis is increased, then thelevel of the hormone hepcidin is decreased. Thus, for instance, lowhepcidin concentrations are observed in patients with absoluteiron-deficiency anemia, or anemias with high erythropoietic activity.See Goodnough 2010. The consequence of a lower hepcidin level is thatmore iron is taken up via the duodenum and more iron is released frommacrophages and from the stores. See Crichton; Fleming. If iron storesare full and/or erythropoiesis is decreased, then more hepcidin issecreted by hepatocytes. Dietary iron absorption, macrophage ironrecycling and release of iron from stores is partly or fully blocked asa result of the binding to the iron exporter ferroportin.

Hepcidin is an acute phase-reactant, and, during inflammation, cytokinesstimulate overproduction of hepcidin. Hence, the feedback is alteredduring inflammation, which is a prominent problem in chronic kidneydisease. As a consequence, the body is not able to sufficiently supplydeveloping red blood cells with iron and erythropoiesis is impaired,although there might be more than enough iron in stores and macrophages,but it is simply not available. This paradoxical situation is calledfunctional iron deficiency, a condition where stored iron is sufficient,but circulating iron is deficient, whereas a situation where the ironcontent of the body is too low (depleted stores) is referred to asabsolute iron deficiency that occurs when iron stores are depleted.

7.1 Iron Model

Assumptions

-   -   25. Five iron pools are considered within the body: a plasma        pool, iron in precursor cells, iron in erythrocytes, iron in the        macrophages and a storage pool.    -   26. Other iron in the body is ignored.    -   27. Hepcidin regulates iron homeostasis.    -   28. Stressed erythropoiesis and depleted stores decrease the        amount of hepcidin released into plasma.    -   29. Large iron stores and inflammation increase hepcidin levels.    -   30. The amount of iron in plasma is proportional to the amount        of transferrin molecules carrying iron. (Free iron is        neglected.)    -   31. There is no differentiation between monoferric and differic        transferrin.    -   32. Iron is only lost via loss of cells, i.e. iron lost via        urine and sweat is neglected.

Iron is a part of hemoglobin which provides the means of O₂ transport tothe tissue. The rate at which iron can be released to plasma fromexisting stores and macrophages may easily limit the rate of irondelivery for hemoglobin synthesis. Therefore, a strong relationshipbetween impaired erythropoiesis and total body iron burden exists. SeeCrichton for an excellent description of iron homeostasis and itseffects.

Iron is a substance which is very tightly controlled from the body asiron overload is toxic. Iron homeostasis can be only achieved by thecontrol of absorption because the human body, unlike other mammals, isnot able to influence the excretion of iron. Under normal circumstancesthe average daily loss of iron is very small (men 1 mg, women 2 mg,because of menstruation).

For the Iron Model, five dynamic iron pools are considered within thebody: a plasma pool, iron in precursor cells, iron in erythrocytes, ironin the macrophages and a storage pool. Hence, ferritin and hemosiderinstores are lumped together. Other iron in the body is ignored. Theamount of iron lost via urine and sweat is neglected. Thus, in thismodel, excretion of iron only takes place when cells are lost. On onehand, this can be due to loss of RBCs, i.e., due to external bleeding,and, on the other hand, there is a daily loss of epithelial cells.Further, one assumes that only the precursor cells in the bone marrowuptake iron, and ineffective erythropoiesis depends on the amount ofiron which can be provided for erythropoiesis. Moreover, the modelaccounts for the loss of iron from erythrocytes during senescence.

The general form of the compartmental Iron Model is as follows:

${{\frac{d}{dt}a_{1}} = {{{k_{41}(H)}a_{4}} + {h_{{iv},1}a_{iv}} + {{k_{{gastro},1}(H)}a_{gastro}} + {{k_{51}(H)}a_{5}} - {k_{15}a_{1}} - {k_{12}a_{1}}}},$${{\frac{d}{dt}a_{2}} = {{k_{12}a_{1}} - f_{24} - f_{23}}},$${{\frac{d}{dt}a_{3}} = {f_{23} - f_{34} - f_{3,{out}}}},$${{\frac{d}{dt}a_{4}} = {f_{24} + {f_{34}a_{3}} - {{k_{41}(H)}a_{4}} - {k_{45}a_{4}}}},$${\frac{d}{dt}a_{5}} = {{k_{15}a_{1}} - {k_{45}a_{4}} - {{k_{51}(H)}{a_{5}.}}}$

Note, in the above equations, the dependence of particular components ont was not explicitly stated to simplify notation. Furthermore, a₁, i=1,. . . , 5, denotes the amount of iron in the compartment i and k_(y), i,j=1,000.5 are the transfer rates from compartment i to compartment j,i.e., k_(y)a_(i) is the rate at which iron is moving from compartment ito compartment j. Further, ƒ_(ij) denote the rates at which iron istransferred from compartment i into compartment j, when it is not of thesimple form k_(ij)α_(i); ƒ_(3,out) is the amount of iron which is lostby bleeding. Additionally, α_(iv) denotes the amount of ironintravenously administered, whereas α_(gastro) is the amount of iron inthe duodenum. Finally, k_(iv,1) and k_(gastro,1)(H) are thecorresponding flow rates. The rate k_(gastro,1) can be described by adecreasing sigmoidal function

${{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + e^{{{- k_{5}}{H(t)}} + c_{5}}} + b_{5}}},$

where a₅, b₅, c₅ and k₅ are positive constants with a₈>b₈. The functionH(t) is the solution of the differential equation (13) below.

Since the time needed to administer iron is very short compared to othertime constants in the system, the following can be taken

k _(iv,1) a _(iv) =a _(iv,total)δ_(t) ₀ .

Here a_(iv,total) is the total amount of iron administered and δ_(t0) isthe Dirac delta function located at t₀.

7.1.1 Hepcidin Feedback

The regulator of iron homeostasis is hepcidin, a peptide hormoneproduced by the liver. Hepcidin negatively regulates the availability ofiron in plasma by direct inhibition of ferroportin. Ferroportin is theiron exporter required for iron egress from iron exporting cells, as forinstance, enterocytes and macrophages. Hepcidin acts as a medium towithhold iron from certain invasive bacteria, as they are not able toproliferate under conditions of insufficient iron supply. The dynamicsof the hepcidin concentration in plasma is calculated according to

$\begin{matrix}{{{\frac{d}{dt}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{H(t)}}}},} & (13)\end{matrix}$

where H_(prod)(t) is the rate at which hepcidin, produced by thehepatocytes, is released into the plasma at time t. TBV is the totalblood volume, and c^(H) _(deg) is the rate at which hepcidin isdegraded. The production rate H_(prod)(t) of hepcidin is obtained as theresult of a feedback loop involving the total precursor cell population,the amount of iron stored and the state of inflammation of the patient,

H _(prod)(t)=ƒ(U(t),a ₅(t),i(t)).

Here

${U(t)} = {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{s\left( {t,\mu^{s},} \right)}d\mu^{s}}}}$

where U(t) describes the total precursor cell population(erythroblasts+reticulocytes). Further, a₅(t) is the amount of ironstored and i(t) is the inflammation status at time t. The feedback ismodeled using a monotonically increasing sigmoidal function

${H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + e^{{{- k_{H}}{C(t)}} + c_{H}}} + b_{H}}$

where a_(H), b_(H),c_(H) and k_(H) are positive constants witha_(H)>b_(H) The function

C(t)=c _(U) U(t)+c _(a) ₅ a ₅(t)+c _(i) i(t)

is an auxiliary equation that increases if stores or inflammation statusrises, and decreases when more RBCs are produced. Here, c_(U), c_(a5)and c_(i) are positive constants.

7.2. Iron Model Linked to the Erythropoiesis Mode.

Assumptions

-   -   33. All types of precursor cells synthesize hemoglobin.    -   34. The rate at which the amount of hemoglobin in precursor        cells changes, depends on the number of transferrin receptors        (TfR) on the cell membrane and the amount of iron in plasma.    -   35. Cells with the same age express the same amount of TfR on        the membrane.    -   36. An erythroblast starts with 300,000 TfR on the cell        membrane.    -   37. This number increases sharply and peaks after 1.5 days at        800,000 TfR.    -   38. After that time, it declines to 100,000 TfR on bone marrow        reticulocytes.    -   39. The number of TfR on bone marrow reticulocytes is constant        (100,000 TfR/cell).    -   40. The rate of apoptosis of precursor cells (erythroblasts AND        reticulocytes) depend on cell age and hemoglobin.    -   41. Iron that enters a precursor cell is incorporated into        hemoglobin, i.e., no free iron in the cell is considered.    -   42. Erythrocytes don't lose iron during senescence.

Some thought has to be given to how the iron model can be linked to thepopulation model. Iron is consumed by precursor cells (erythroblasts andbone marrow reticulocytes) and used for hemoglobin synthesis. The amountof iron that is taken up by a cell depends on the number of transferrinreceptors (TfR) expressed on the membrane and the concentration of ironin plasma. In order to simplify notation, the minimal cell age ofprecursor cells p min is set to zero, i.e., the cell age forerythroblasts is shifted from 8=μ^(r) _(min)≤μ^(r)≤μ^(r) _(max)=13 to0=μ^(r) _(min)≤μ^(r)≤μ^(r) _(max)=5 and the cell age for bone marrowreticulocytes consequently is shifted from 13=μ^(s) _(min)≤μ^(s)≤μ^(s)_(max)=15 to 5=μ^(s) _(min)≤μ^(s)≤μ^(s) _(max)=7.

The contact rate between the TfRs of a cell and transferrin molecules inplasma carrying iron is governed by the law of mass action and thus isgiven by:

$\begin{matrix}{{d_{1}{\tau(t)}\frac{a_{1}(t)}{TBV}},} & (14)\end{matrix}$

where d₁ is a rate constant that differs for different patients, as itdepends on the volume of the active bone marrow. Further, τ(t) is thecumulative amount of TfR on cell membranes of precursor cells(erythroblasts and bone marrow reticulocytes) at time t and α₁(t) is theconcentration of iron in plasma at time t. Note that this is at the sametime the rate with which iron leaves the plasma compartment and entersthe precursor cells. The cumulative sum of TfR expressed on allprecursor cells τ(t) is determined by

$\begin{matrix}{{{\tau(t)} = {{\int\limits_{\mu_{\min}^{T}}^{\mu_{\max}^{T}}{{\varphi\left( \mu^{r} \right)}{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{\varphi\left( \mu^{s} \right)}{s\left( {t,\mu^{s},} \right)}d\mu^{s}}}}},} & (15)\end{matrix}$

where φ(μ) is the function describing the number of TfR on a cell withcell age μ, r(t,μ^(r)) is the density of erythroblasts, and s(t,μ^(s))is the density of bone marrow reticulocytes.

7.2.1. Function for the Transferrin Receptors.

Assumptions 35-39 govern the development of the transferrin receptors.One ensures that Assumptions 35-38 are met by choosing a polynomial ofdegree 3 φ(μ)=a₄u³+b₄μ²+c₄μ+d₄

and choosing the coefficients a₄, b₄, c and d₄ such that the function φsatisfies the following conditions:

φ(0)=30,φ(1.5)=80,φ(1.5)=0,φ(5)=10.  (16)

For computational reasons, the number of TfR is multiplied by 10⁻⁴. Thenthe first condition in eq. (16) coincides with Assumption 36. The secondand third condition reflect Assumption 37 and, with the last condition,one ensures that there is a smooth transition in the number of TfR fromerythroblastic cells to reticulocytes. Deriving and solving the linearsystem that arises from eq. (16) gives a₄=3.3016 b₄=−32.127 d₄=30.Further, using Assumption 38, the function ø can be extended to bonemarrow reticulocytes. Altogether, one has

$\begin{matrix}{{\varphi(\mu)} = \left\{ \begin{matrix}{{{3.3016\mu^{3}} - {32.127\mu^{2}} + {74.0952\mu} + 30},{{for}{{\mu\epsilon}\left\lbrack {\mu_{\min}^{r},\mu_{\max}^{r}} \right\rbrack}},} \\{10,{{for}{{\mu\epsilon}\left\lbrack {\mu_{\min}^{s},\mu_{\max}^{s}} \right\rbrack}}}\end{matrix} \right.} & (17)\end{matrix}$

7.2.2. Accumulation of Hemoglobin

Red blood cells carry hundreds of hemoglobin molecules which enablesthem to carry oxygen. Hemoglobin has four subunits each containing aheme group. An iron atom is embedded in the middle of the heme group.Hence, in one hemoglobin molecule four iron atoms are incorporated. Inorder to keep track of the accumulation of hemoglobin molecules in theprecursor cells, the development of the TfR and the amount of ironavailable at certain times have to be considered. A consequence ofAssumption 35 is that cells with the same cell age all carry the sameamount of hemoglobin. The amount of hemoglobin a cell with ū accumulatedcan vary over time t but is invariant among cells with age ū at aspecific time t. First, one focuses on the erythroblast populationclass, where one has the special case that the maturation velocity v≡1,i.e., the cell age and the time for how long a cell has been anerythroblast, and thus had time to accumulate hemoglobin, coincide.Hereinafter, the actual age of the cell, i.e. the time that elapsedsince the cell became a precursor cell, will be called calendric age.This is to emphasize the difference between the cell age and the actualage, i.e. the calendric age.

The rate with which hemoglobin increases in a cell is given by

(18)

${{\frac{d}{dt}\left( {h\left( {t,{\mu(t)}} \right)} \right)} = {{{\frac{\partial}{\partial t}{h\left( {t,{\mu(t)}} \right)}} + {{\frac{\partial}{\partial\mu}{h\left( {t,{\mu(t)}} \right)}}{v\left( {E(t)} \right)}}} = {d{\varphi\left( {\mu(t)} \right)}\frac{a_{1}(t)}{TBV}}}},$

where h(t,μ(t)) denotes the amount of hemoglobin in a cell with μ (t) attime t. Note that, in most cases, the dependence of the cell age μ ontime t is not explicitly stated, because there is no necessity to do so.Here, one has to consider this fact, as it is important for thederivation of the function h. Further, d is a rate constant which isdifferent from d₁ because here one considers hemoglobin instead of iron,see eq. (14). The ratio between d and d₁ is given by the fact that eachhemoglobin molecule consists of four iron atoms. The function φ(μ) isdetermined by eq. (17), α₁(t) describes the amount of iron in the plasmacompartment and TBV denotes the total blood volume. For a detaileddescription of how the amount of hemoglobin in a cell with age μ at timet can be calculated from eq. (18).Maturation velocity v≡1

Cells in the erythroblast population class in the Erythropoiesis Modelsatisfy the condition that cell age and calendric age of the cellcoincide, see Section 6. This makes it easier to keep track of theaccumulation of hemoglobin. In order to calculate the amount ofhemoglobin a cell with cell age μ at time t carries, one has to considerthat the cell entered the population class at time t−μ+μ_(min), whereμ_(min) denotes the minimal cell age of the population class. Further,the amount of TfR on the cell membrane of the cell that entered thepopulation at time t−μ+μ_(min) was φ(μ_(min)) and the amount of iron inplasma at that time was α₁(t−μ+μ_(min)). With these considerations, theamount of hemoglobin a cell with cell age .mu. at time t carries isgiven by the following integral

${h\left( {\overset{\_}{t},\overset{\_}{\mu}} \right)} = {\frac{d}{TBV}{\int\limits_{\overset{\_}{t} - \overset{\_}{\mu} + \mu_{\min}}^{\overset{\_}{t}}{{\varphi\left( {s - \overset{\_}{t} + \overset{\_}{\mu}} \right)}{a_{1}(s)}{ds}}}}$

Note that the calculation makes use of the history of the availableamount of iron in plasma from t−μ+μ_(min) to t.

For the erythroblast population (the cell age was transformed such thatμ^(r) _(min)≡0) the above equation simplifies to

$\begin{matrix}{{h\left( {t,\mu^{r}} \right)} = {\frac{d}{TBV}{\int\limits_{t - \mu^{r}}^{t}{{\varphi\left( {s - t + \mu^{r}} \right)}{a_{1}(s)}{{ds}.}}}}} & (19)\end{matrix}$

Maturation velocity v≡v(E(t))

In the case that cells have a variable maturation velocity, it gets moredifficult to keep track of the accumulation of hemoglobin. First, onehas to calculate at which time φ {circumflex over (t)} a cell that hascell age μ time t entered the population class. Therefore, the followingequation has to be solved with regard to {circumflex over (t)}

$\overset{\_}{\mu} = {\mu_{\min} + {\int\limits_{\hat{t}}^{\overset{\_}{t}}{{v\left( {E(s)} \right)}{ds}}}}$

i.e., one has to (numerically) find the root of the following function

$\begin{matrix}{{f\left( {\overset{\_}{t},\overset{\_}{\mu}} \right)} = {\mu_{\min} - \overset{\_}{\mu} + {\int\limits_{\hat{t}}^{\overset{\_}{t}}{{v\left( {E(s)} \right)}{{ds}.}}}}} & (20)\end{matrix}$

The value η=t−{circumflex over (t)} reflects the calendric age of thecell. (It represents the time that passed since the cell entered thepopulation class and does not necessarily coincide with the cellage/maturity of the cell.)

In the Erythropoiesis model presented in Section 6, only thereticulocyte population has a variable maturation velocity. Assumption39 states that cells in this specific population class express aconstant number of TfR on their membranes. For the sake of simplicity,the special case φ≡const is considered here, and the general case whenφ=(μ) is not considered.

Thus, the amount of hemoglobin in a reticulocyte with calendric ageη^(s)ϵ[5,8] at time t is given

$\begin{matrix}{{h\left( {t,\eta^{s}} \right)} = {{h\left( {{t - \eta^{s} + \mu_{\min}^{s}},\mu_{\max}^{r}} \right)} + {\varphi\frac{d}{TBV}{\int\limits_{t - \eta^{s} + \mu_{\min}^{s}}^{t}{{a_{1}(s)}{{ds}.}}}}}} & (21)\end{matrix}$

The calendric age of the cell has been derived using the solution of eq.(20) and h(t−η^(s)+μ^(s) _(min), μ^(r) _(max)) is derived from eq. (19).

7.2.3 Rate of Apoptosis of Precursor Cells

The construction of the rate of apoptosis of precursor cells is inspiredby the fact that cells that do not synthesize hemoglobin sufficientlyare more likely to die. See Schmidt.

Hereinafter, one is going to refer to cell age of the erythroblasts andthe calendric age (see eq. 20) of the reticulocytes, respectively, as η.Consequently, ηϵ[0,8] and denotes the time that elapsed since a cellbecame a precursor cell.

First, a “reference hemoglobin accumulation curve” denoted by h*(η) iscreated using eqs. (19) and (21). The construction is based on anaverage male healthy adult, assuming that: the person's erythropoiesisand iron homeostasis is in steady state, the total blood volume is 5liters, and the amount of iron in plasma is 3 mg, i.e., α₁≡3.

Further, healthy persons produce normocytic red blood cells and theamount of Hgb per cell lies within 26-32 pg. The contact rate d ischosen such that, under the above assumptions and by applying eqs. (19)and (21), respectively, the hemoglobin content of a cell lies wellwithin this range when the cell leaves the reticulocyte population andenters the blood stream. Note that the label ‘days since cell became aprecursor cell’ is similar to the calendric age .eta. and forerythroblast (day 0-5) coincides with the cell age, whereas for thereticulocytes the cell age and the calendric age might differ.Reticulocytes may stay in the bone marrow 1-3 days, but they all leavethe bone marrow when they reach maximal maturity, i.e., when their cellage turns 2. This means that reticulocytes might have differenttimespans to gather iron and synthesize hemoglobin—depending on thematuration velocity v(E(t)), see Section 6.

Next, the reference curve is shifted, by slightly reducing orincreasing, respectively, the amount of iron in plasma (α₁≡2.75 mg andα₁≡3.25 mg, respectively). The consequence being that a corridordevelops around the reference curve. The next step is to designfunctions in a way that cells that leave this corridor are more likelyto die, the farther they diverge from it, the higher the rate ofapoptosis gets. A family of functions is created that describes thisincrease in mortality as cells diverge from the reference curve

$\begin{matrix}{{\aleph\left( {h,\eta} \right)} = {\left( \frac{c_{\aleph}}{\left( {\eta + 1} \right)^{3}} \right)\left( {h - {h^{*}(\eta)}} \right)^{4}}} & (22)\end{matrix}$

where

is a constant, h denotes the amount of hemoglobin that a cell carries atcell age/calendric-age η and h* is the reference hemoglobin value at theage η. Polynomial functions of degree 4 are used to design thisfunction, because one of their characteristics is that they have a flatbottom (i.e. apoptosis within the corridor is low) and then increasesharply. Finally, to get the rate of apoptosis α for the precursorcells, a discrete set of N+1 equidistant data points μ₁, i=0, N aredefined along μ∈[μ_(min) ^(r),μ_(max) ^(s)] at which the followingcalculations are done calculate the amount of Hgb h_(i) in cells of ageμ_(i) using eq. (19), (20) and (21), calculate the values of thefunction

at the grid points η_(i) and the corresponding Hgb value h_(i)interpolate the resulting values using cubic splines. The result is thefunction α(h,η) describing the rate of apoptosis of precursor cells.Note that the above calculations have to be repeated at every time step.

7.3. Detailed Description of the Iron Compartments. 7.3.1. Plasma(Compartment 1).

The rate at which iron enters erythroid cells at time t is given by

d ₁ a ₁(t)τ(t)

where the number of transferrin receptors τ on all precursor cells is

${\tau(t)} = {d_{1}\left( {{\int\limits_{\mu_{\min}^{r}}^{\mu_{\max}^{r}}{{\varphi\left( \mu^{r} \right)}{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\int\limits_{\mu_{\min}^{s}}^{\mu_{\max}^{s}}{{\varphi\left( \mu^{s} \right)}{s\left( {t,\mu^{s}} \right)}d\mu^{s}}}} \right)}$

where φ(μ) is given by eq. (17). Further, r(t,μ^(r)) and s(t,μ^(s))denote the population densities of hemoglobin synthesizing cells,erythroblasts and reticulocytes, respectively. At the same time this isthe rate at which iron leaves the plasma compartment and enters theprecursor cell compartment.

Only a very small amount of iron can be found in this pool (about 3 mg;.ltoreq.1%) bound to transferrin, but it is a very important compartmentbecause of its rapid turnover rate. See Finch 1982. Iron normally turnsover at least 10 times each day. See Williams Hematology. Inflow intothe plasma compartment comes from macrophages, the storage pool, ironabsorbed through the duodenum, and via intravenously administered iron.Outflow leaves the plasma compartment to the precursor cells and thestorage compartment. The flow rates k₄₁,k_(gastro,1) and k₅₁ aredependent on the hepcidin concentration H(t) at time t. They decreasewhen the hepcidin level rises and vice versa. Inflammation and full ironstores increase the release of hepcidin from the hepatocytes in theliver, whereas a stressed erythropoiesis and depleted stores suppressthe production of hepcidin. The outflow to the precursor cellcompartment k₁₂α₁ is altered by the amount of iron needed from the cellsto synthesize hemoglobin. As erythroid cells are only able to take upiron through transferrin receptors, their need for iron can beassociated with the number of receptors expressed on precursor cells.Thus, the rate k₁₂ is given by

$k_{12} = {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\varphi\left( \mu^{r} \right){r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\mu_{\min}^{s}}{\int\limits^{\mu_{\max}^{s}}}{\varphi\left( \mu^{s} \right){s\left( {t,\mu^{s},} \right)}d\mu^{s}}}}$

For an explanation of the various terms see the statements abouterythroblasts in Section 6.2. Hence, the change of amount of iron in theplasma compartment is

${\frac{d}{dt}a_{1}} = {{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} + {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{1}\left( {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\varphi\left( \mu^{r} \right){r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\mu_{\min}^{s}}{\int\limits^{\mu_{\max}^{s}}}{\varphi\left( \mu^{s} \right){s\left( {t,\mu^{s},} \right)}d\mu^{s}}}} \right)}{{a_{1}(t)}.}}}$

The transfer rate k₄₁(H(t)) and k₅₁(H(t)) can be described bymonotonically decreasing sigmoidal functions

$\begin{matrix}{{k_{41}\left( {H(t)} \right)} = {\frac{a_{6} - b_{6}}{1 + e^{{{- k_{6}}{H(t)}} + c_{6}}} + b_{6,}}} & (23)\end{matrix}$ and $\begin{matrix}{{k_{51}\left( {H(t)} \right)} = {\frac{a_{7} - b_{7}}{1 + e^{{{- k_{7}}{H(t)}} + c_{7}}} + b_{7,}}} & (24)\end{matrix}$

respectively. The parameters a₆, b₆, c₆, k₆, a₇, b₇, c₇ and k₇ arepositive constants with a₆>b₆ and a₇>b₇. The hepcidin level H(t) at timet is determined by eq. (13).

7.3.2. Precursor Cells (Compartment 2).

Inflow into the precursor cell compartment comes from plasma. Oneoutflow of this compartment leaves to the erythrocytes compartment anddescribes the amount of hemoglobin carried from reticulocytes, whichenter the blood stream. Thus, the flow ƒ₂₃(⋅) from compartment 2 tocompartment 3, is defined by the cells that leave the bone marrow andthe amount of hemoglobin these cells are carrying. First, one has tocompute how long cells that leave the bone marrow at time t haveactually stayed there as precursor cells. Hence, one needs to solveƒ({circumflex over (t)}, μ_(max) ^(s)) and one denotes the root of eq.(20) by η^(s) _(max). Note that η^(smax) might change over time asreticulocytes are released in the blood stream faster or more slowly.

k ₂₃ a ₂=ƒ₂₃(v ^(s)(E(t))s(t,μ _(max) ^(s)),h(t,μ _(max) ^(s))),

where v^(s)(E(t))s(t,μ^(s) _(max)) determines how many reticulocytes areleaving the bone marrow at a certain time t and the functionh^(s)(t,η^(s) _(max)) indicates how much hemoglobin they carry.v^(s)(E(t)) is the maturation velocity of reticulocytes depending on EPOand s(t,μ^(s) _(max)) is the density of the reticulocyte population withmaximal maturity μ^(s) _(max) at time t. The function h^(s)(t,μ^(s)_(max)) is given by eq. (21).

The second outflow of this compartment leaves to the macrophages and isdefined by the precursor cells that die. Hence, the flow fromcompartment 2 to compartment 4 is defined by the amount of iron thatdying cells carry, i.e.

$\begin{matrix}{{k_{24}a_{4}} = {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\alpha\left( {h,\mu^{r}} \right){h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\eta_{\min}^{s}}{\int\limits^{\eta_{\max^{(t)}}^{s}}}{\alpha\left( {h,\eta^{s}} \right){h\left( {t,\eta^{s}} \right)}{s{}\left( {t,{\mu_{\min} + {\underset{t - \eta^{s}}{\int\limits^{t}}{{v\left( {E(s)} \right)}{ds}}}}} \right)}d{\eta^{s}.}}}}} & (25)\end{matrix}$

For the first part of the above equation, η and μ^(r) are exchangeablewithout any difficulties as, for erythroblasts, both describe the cellage of the cell, see Section 6.2.

Altogether the change of amount of iron in the iron-precursorscompartment is given by

${\frac{d}{dt}{a_{2}(t)}} = {{{d_{0}\left( {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\varphi\left( \mu^{r} \right){r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\mu_{\min}^{s}}{\int\limits^{\mu_{\max}^{s}}}{\varphi\left( \mu^{s} \right){s\left( {t,\mu^{s},} \right)}d\mu^{s}}}} \right)}{a_{1}(t)}} - {f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s}} \right)}},{h^{s}\left( {t,\mu_{\max}^{s}} \right)}} \right)} - {\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{{\alpha\left( {h,\mu^{r}} \right)}{h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} - {\underset{\eta_{\min}^{s}}{\int\limits^{\eta_{\max^{(t)}}^{s}}}{\alpha\left( {h,\eta^{s}} \right){h\left( {t,\eta^{s}} \right)}s\left( {t,{\mu_{\min} + {\underset{t - \eta^{s}}{\int\limits^{t}}{{v\left( {E(s)} \right)}{ds}}}}} \right)d{\eta^{s}.}}}}$

Note that it is not necessary to determine the function ƒ₂₃ as, in fact,there is no need to explicitly observe the rate of change of theprecursor cell compartment. The amount of iron α²(t) in this compartmentat time t can be directly computed from eqs. (19), (20), (21), and thepopulation equations for the erythroblasts and reticulocytes. It isdefined by the total number of precursor cells and the amount ofhemoglobin they carry:

${a_{2}(t)} = {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\alpha\left( {h,\mu^{r}} \right){r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\eta_{\min}^{s}}{\int\limits^{\eta_{\max^{(t)}}^{s}}}{\alpha\left( {h,\eta^{s}} \right){h\left( {t,\eta^{s}} \right)}{s{}\left( {t,{\mu_{\min} + {\underset{t - \eta^{s}}{\int\limits^{t}}{{v\left( {E(s)} \right)}{ds}}}}} \right)}}}}$

However, it is important to know the flux of hemoglobin to themacrophage compartment, which arises from dying cells, because this flux(see eq. (25)) is a source term in the macrophage compartment.

7.3.3. Erythrocytes (Compartment 3)

Inflow into the erythrocytes compartment comes from the precursor cellcompartment, and iron leaves the compartment to the macrophagescompartment (phagocytosis of senescent cells, random break-down of cellsand neocytolysis)

$\begin{matrix}{{{k_{34}a_{3}} = {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}d\mu^{m}}}},} & (26)\end{matrix}$

or is excreted via bleeding

${k_{3,{out}}a_{3}} = {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha_{bleed}^{m}(t)}{m\left( {t,\mu^{m}} \right)}d\mu^{m}}}$

An erythrocyte of age μ^(m) left the reticulocyte population at timet−p^(m) and, making use of Assumption 42, the amount of hemoglobin itcarries can be computed using eq. (21). Further, α^(m) denotes themortality rate of erythrocytes due to phagocytosis of senescent cells,random break-down of cells, and neocytolysis. Furthermore, h is theamount of hemoglobin carried by a cell, m(t, μ^(m)) is the populationdensity of erythrocytes with cell age μ^(m) at time t, and α^(m)_(bleed) is the rate with which RBCs are lost due to bleeding.

Thus, the change of the amount of iron in the erythrocyte compartment is

${\frac{d}{dt}{a_{3}(t)}} = {{f_{23}\left( {{{v^{s}\left( {E(t)} \right)}{s\left( {t,\mu_{\max}^{s}} \right)}},{h^{s}\left( {t,\mu_{\max}^{s}} \right)}} \right)} - {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}d\mu^{m}}} - {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha_{bleed}^{m}(t)}{m\left( {t,\mu^{m}} \right)}d{\mu^{m}.}}}}$

Similar to compartment 2, it is not necessary to observe the rate ofchange of the erythrocyte cell compartment to compute the amount of ironα₃(t) at time t. The quantity α₃(t) can be directly computed from thepopulation class and is defined by the total number of circulating redblood cells and the amount of hemoglobin they carry.

${a_{3}(t)} = {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{m\left( {t,\mu^{m}} \right)}d\mu^{m}}}$

However, the term k₃₄α₃(t) needs to be known explicitly (see eq. (26)).This flux, which is the rate of hemoglobin released by erythrocytes,appears as a source term in the macrophage compartment.

7.3.4 Macrophages (Compartment 4)

Inflow into the macrophages compartment comes from the erythrocytescompartment (erythrocytes phagocytosed by macrophages, iron set free byrandom break-down of RBCs and shedding of hemoglobin containingvesicles) and the precursor cells compartment. Outflow of thiscompartment leaves to the plasma and to the storage compartment.Therefore, the change of iron in the compartment is

${{\frac{d}{dt}{a_{4}(t)}} = {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\alpha\left( {h,\mu^{r}} \right){h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\eta_{\min}^{s}}{\int\limits^{\eta_{\max^{(t)}}}}{\alpha\left( {h,\eta^{s}} \right){h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\underset{t - \eta^{s}}{\int\limits^{t}}{{v\left( {E(s)} \right)}{ds}}}}} \right)}d\eta^{s}}} + {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}d\mu^{m}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}},$

where k₄₁(H(t)) is given by eq. (23).

7.3.5. Storage (Compartment 5)

Inflow into the storage compartment comes from the plasma and themacrophages compartment and iron leaves this compartment to the plasmacompartment. Hence, the change of iron in the storage compartment is

${\frac{d}{dt}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} + {k_{45}{a_{4}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}$

where k₅₁(H(i)) is given in eq. (24).7.4. Mathematical Model lion Model

For the convenience of the reader, all equations of the Iron Model arecollected in this section:

${{\frac{d}{dt}{a_{1}(t)}} = {{{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} + {{a_{{iv},{total}}(t)}{\delta_{t_{0}}(t)}} + {{k_{{gastro},1}\left( {H(t)} \right)}{a_{gastro}(t)}} + {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}} - {k_{15}{a_{1}(t)}} - {{d_{1}\left( {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\varphi\left( \mu^{r} \right){r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\mu_{\min}^{s}}{\int\limits^{\mu_{\max}^{s}}}{\varphi\left( \mu^{s} \right){s\left( {t,\mu^{s},} \right)}d\mu^{s}}}} \right)}{a_{1}(t)}}}},$${a_{2}(t)} = {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\alpha\left( {h,\mu^{r}} \right){r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\mu_{\min}^{s}}{\int\limits^{\eta_{\max^{(t)}}^{s}}}{\alpha\left( {h,\eta^{s}} \right){h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\underset{t - \eta^{s}}{\int\limits^{t}}{{v\left( {E(s)} \right)}{ds}}}}} \right)}d\eta^{s}}}}$${a_{3}(t)} = {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{m\left( {t,\mu^{m}} \right)}d\mu^{m}}}$${{\frac{d}{dt}{a_{4}(t)}} = {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{\alpha\left( {h,\mu^{r}} \right){h\left( {t,\mu^{r}} \right)}{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\mu_{\min}^{s}}{\int\limits^{\eta_{\max^{(t)}}^{s}}}{\alpha\left( {h,\eta^{s}} \right){h\left( {t,\eta^{s}} \right)}{s\left( {t,{\mu_{\min} + {\underset{t - \eta^{s}}{\int\limits^{t}}{{v\left( {E(s)} \right)}{ds}}}}} \right)}d\eta^{s}}} + {\underset{\mu_{\min}^{m}}{\int\limits^{\mu_{\max}^{m}}}{{h\left( {{t - \mu^{m}},\eta_{\max}^{s}} \right)}{\alpha^{m}\left( {{E(t)},\mu^{m}} \right)}{m\left( {t,\mu^{m}} \right)}d\mu^{m}}} - {{k_{41}\left( {H(t)} \right)}{a_{4}(t)}} - {k_{45}{a_{4}(t)}}}},$${{\frac{d}{dt}{a_{5}(t)}} = {{k_{15}{a_{1}(t)}} - {{k_{51}\left( {H(t)} \right)}{a_{5}(t)}}}},$${\frac{d}{dt}{H(t)}} = {{\frac{1}{TBV}{H_{prod}(t)}} - {c_{\deg}^{H}{H(t)}}}$

where the initial values are given by

a ₁(0)=a _(1,0),

a ₂(0)=a _(2,0),

a ₃(0)=a _(3,0),

a ₄(0)=a _(4,0),

a ₅(0)=a _(5,0),

H(0)=H ₀.

Moreover, a number of auxiliary equations are needed

${{k_{{gastro},1}\left( {H(t)} \right)} = {\frac{a_{5} - b_{5}}{1 + e^{{{- k_{5}}{H(t)}} + c_{5}}} + b_{5}}},$${{k_{41}\left( {H(t)} \right)} = {\frac{a_{6} - b_{6}}{1 + e^{{{- k_{6}}{H(t)}} + c_{6}}} + b_{6}}},$${{k_{51}\left( {H(t)} \right)} = {\frac{a_{7} - b_{7}}{1 + e^{{{- k_{7}}{H(t)}} + c_{7}}} + b_{7}}},$

where the hepcidin feedback is given by

${{H_{prod}(t)} = {\frac{a_{H} - b_{H}}{1 + e^{{{- k_{H}}{C(t)}} + c_{H}}} + b_{H}}},$C(t) = −c_(U)U(t) + c_(a₅)a₅(t) + c_(i)i(t),${U(t)} = {{\underset{\mu_{\min}^{r}}{\int\limits^{\mu_{\max}^{r}}}{{r\left( {t,\mu^{r}} \right)}d\mu^{r}}} + {\underset{\mu_{\min}^{s}}{\int\limits^{\mu_{\max}^{s}}}{{s\left( {t,\mu^{s}} \right)}d\mu^{s}}}}$

Furthermore, we have

${\varphi(\mu)} = \left\{ \begin{matrix}{{{3.3016\mu^{3}} - {32.127\mu^{2}} + {74.0952\mu} + 30},{{{for}\mu} \in \left\lbrack {\mu_{\min}^{r},\mu_{\max}^{r}} \right\rbrack}} \\{10,{{{for}\mu} \in \left\lbrack {\mu_{\min}^{s},\mu_{\max}^{s}} \right\rbrack}}\end{matrix} \right.$${h\left( {t,\mu^{r}} \right)} = {\frac{d}{TBV}{\overset{t}{\int\limits_{\iota - \mu^{r}}}{{\varphi\left( {s - t + \mu^{r}} \right)}{a_{1}(s)}{ds}}}}$${h\left( {t,\eta^{r}} \right)} = {{h\left( {{t - \eta^{s} + \mu_{\min}^{s}},\mu_{\max}^{s}} \right)} + {\phi\frac{d}{TBV}{\overset{t}{\int\limits_{t - \eta^{s} + \mu_{\min}^{s}}}{{a_{1}(s)}{{ds}.}}}}}$

The calendric age η of cells is calculated using

${{f\left( {\overset{\_}{t},\overset{\_}{\mu}} \right)} = {\mu_{\min} - \overset{\_}{\mu} + {\overset{\overset{\_}{t}}{\int\limits_{\hat{t}}}{{v\left( {E(s)} \right)}{ds}}}}},$

and η=t−{circumflex over (t)}.

Further, to create the function α(hη) (for η either μ^(r) or μ^(s) canbe inserted) a discrete set of N+1 equidistant data points μ_(i), . . ., N, along μ∈[μ_(min) ^(r),μ_(max) ^(s)] is defined and at every timestep the following calculations are done

calculate the amount of Hgb h_(i) in cells of age μ_(i) using theformula for η, ƒ(t,μ), h(t,μ^(r)) and h(t,μ^(s)),

calculate the values of the function

at the grid points η_(i) and the corresponding Hgb value h_(i)

${\aleph\left( {h,\eta} \right)} = {\left( \frac{c_{\aleph}}{\left( {\eta + 1} \right)^{3}} \right)\left( {h - {h^{*}(\eta)}} \right)^{4}}$

interpolate the resulting values using cubic splines.

Derivation of Population Equations

A.1. Population Equations with One Structural Variable.

Let u(t,⋅) denote the density of the population at time t, i.e., for0≤a<b is the number of cells with maturity level ∫_(a) ^(b)u(t, x)dx.The rate of change is

(A.1)

${{\frac{d}{dt}{\int_{a}^{b}{{u\left( {t,x} \right)}{dx}}}} = {\int_{a}^{\overset{˙}{o}}{{u_{1}\left( {t,x} \right)}{dx}}}},{t \geq 0},{0 \leq a < {b.}}$

Here, u is assumed to be continuously differentiable with respect to t.The rate of change with respect tot of the number of cells with maturitylevel in [a,b] is the sum of several components: the flux of cellsacross the boundary points a and h at time t because of changingmaturity; the rate at which cells with maturity in [a,b] die at time t;the rate at which cells with maturity in [a,b] proliferate.

a) Flux Across a and b.

The flux of cells across maturity level is denoted x at time t by q(t,x)and assume that the maturity of a cell is changing according to

{dot over (x)}(t)=v(t,x(t)),t≥0

where the “velocity” v is continuous in t and Lipschitzean with respectto x, t≥0, 0≥0. Moreover, it is assumed that v(t,x)>0, i.e., thematurity of a cell is always strictly increasing.

Then the flux of cells across a maturity level x is given by

q(t,x)=v(t,x)u(t,x),t≥0,x≥0

Consequently, the net flux of cells at time t into the interval [a,b] is

(A.2)

${{q\left( {t,a} \right)} - {q\left( {t,b} \right)}} = {{{{v\left( {t,a} \right)}{u\left( {t,a} \right)}} - {{v\left( {t,b} \right)}{{\iota\iota}\left( {t,b} \right)}}} = {- {\int_{a}^{\overset{˙}{o}}\left( {{\left( {{v\left( {t,x} \right)}{u\left( {t,x} \right)}} \right)_{x}{dx}},{t \geq 0.}} \right.}}}$

b) Mortality.

Let x_(i) ^(N)=a+Δx, i=0, . . . , N with Δx=(b−a)/N be a uniform mesh on[a,b]. By the mean value theorem for integrals, the number of cells withmaturity level xϵ[x_(i) ^(N),x_(i+1) ^(N)] at time t for somex_(i)*ϵ[x_(i) ^(N),x_(i+1) ^(N) is given by

u(t,x _(i)*)Δx.

The rate at which cells with maturity in [x^(N) _(i),x^(N) _(i+1)] dieat time t is assumed to be proportional to the number of cells withmaturity in this interval and therefore is given by

α(t,x _(i)**)u(t,x _(i)*)δx,

where α(t,x)≥0, t>0, x>0, is the mortality rate which is assumed to be acontinuous function and xi**ϵ[_(xNi),x^(N) _(i)+1]. Summing i=0, . . . ,N−1 and taking N→∞ yields

$\begin{matrix}{{\int_{a}^{b}{{\alpha\left( {t,x} \right)}{u\left( {t,x} \right)}{dx}}},{t \geq 0},} & \left( {A\text{.3}} \right)\end{matrix}$

for the rate at which cells with maturity level in [a,b] die at time t.

c) Proliferation.

The assumption concerning proliferation is that a cell with maturitylevel x dividing at time t gives rise to two daughter cells withmaturity x. By analogous considerations as in the case of mortality, thefollowing expression is obtained for the rate at which cells enter theinterval [a,b] at time t because of proliferation:

$\begin{matrix}{{{{2{\int_{a}^{b}{{\beta\left( {t,x} \right)}{u\left( {t,x} \right)}{dx}}}} - {\int_{a}^{b}{{\beta\left( {t,x} \right)}{u\left( {t,x} \right)}d}}} = {\int_{a}^{b}{{\beta\left( {t,x} \right)}{u\left( {t,x} \right)}{dx}}}},{t \geq {0.}}} & \left( {A\text{.4}} \right)\end{matrix}$

where β(t,x)≥0 t≥0, x≥0, is the proliferation rate in the populationwhich is assumed to be continuous.

d) The Model.

Using (A.1)-(A.4) yields

∫_(a)^(b)u_(t)(t, x)dx = −∫_(a)^(b)((v(t, x)ιι(t, x))_(x)dx − ∫_(a)^(b)α(t, x)u(t, x)dx +  ∫_(a)^(b)β(t, x)u(t, x)dx, t ≥ 0, x ≥ 0.

Assuming that the integrands in the above integrals are continuous andobserving that the interval [a, b]⊂R₊ is arbitrary yields

(A.5)

u _(t)(t,x)+((v(t,x)u(t,x))_(x)=(β(t,x)−α(t,x))u(t,x),t≥0,x≥0.

In addition to this equation, an initial condition needs to beprescribed as

u(0,x)=u ₀(x),x≥0

and a boundary condition for x=0. Assuming v(t,0)>0, t≥0, this boundarycondition is written as a flux condition at x=0:

v(t,0)u(t,0)=ƒ₀(t),t>0.

A.2. Population Equations with Several Structural Variables.

Models for populations structured by several attributes can be derivedanalogously to those for populations with just one attribute. Forsimplicity of notation, populations are restricted to three structuralvariables and cell populations with the attributes cell age x, amount ofhemoglobin in the cell y, and amount of transferrin receptors on thecell membrane z are considered.

Let u(t,x,y,z) be the density of the cell population and Ω be a subsetof the positive cone R³ ₊ in R³ which is the admissible set for theattributes. Then the number of cells with attributes (x,y,z)ϵΩ at time tis given by

∫_(Ω)u(t, x, y, z)dω.

Of interest is the rate at which the number of cells with attributes inΩ at time t changes, i.e., in the quantity

${{\frac{d}{dt}{\int_{\Omega}{{u\left( {t,x,y,z} \right)}d\omega}}} = {\int_{\Omega}{{u_{1}\left( {t,x,y,z} \right)}d\omega}}},$

where it is assumed that u is continuously differentiable with respectto t. The rate at which the number of cells with attributes in Ω changesat time t is composed of several components: the flux of cells across ∂Ωbecause the attributes change; the rate at which cells with attributesin Ω die at time t; the rate at which cells enter or leave Ω at time tbecause of proliferation.

a) Flux Across ∂Ω.

The attributes of a cells change according to

{dot over (x)}(t)=ƒ(t,x(t),y(t),z(t)),

{dot over (y)}(t)=g(t,x(t),y(t),z(t)),

ż(t)=h(t,x(t),y(t),z(t)),  (A.7)

where ƒ, g, h are functions R₊xR³ ₊→R which are continuous with respectto t and Lipschitzean with respect to x, y, z. Then the flux of cells attime t and at (x, y, z)∈R₊ ³ is given by

${F\left( {t,x,y,z} \right)} = {{u\left( {t,x,y,z} \right)}\begin{pmatrix}\begin{matrix}{f\left( {t,x,y,z} \right)} \\{g\left( {t,x,y,z} \right)}\end{matrix} \\{h\left( {t,x,y,z} \right)}\end{pmatrix}}$

Consequently the flux of cells across ∂Ω at time t is given by

∫_(∂Ω)

F(t,x,y,z),v(x,y,z)

dσ,  (A.8)

where v(x,y,z) denotes the unit outside normal to an at the point(x,y,z). By Gauss' integral theorem we get

(A.9)

∫_(∂Ω)⟨F((t, x, y, z), v(x, y, z)⟩dσ = ∫_(Ω)divF(t, x, y, z)dω = ∫_(Ω)((f(t, x, y, z)u(t, x, y, z))_(x) + (g(t, x, y, z)u(t, x, y, z))_(y) + (h(t, x, y, z)u(t, x, y, z))_(z))dω.

b) Mortality.

Let Ω=∪_(i−1) ^(N) Ω_(i) such that diam Ω_(i)<ϵ, i=1, . . . , N. By themean value theorem for integrals, the number of cells with attributes inΩ_(i) at time t is given by μ(t, x_(i)*, y_(i)*z_(i)*)ΔΩ_(i), whereΔΩ_(i) denotes the volume of Ω_(i) and (x_(i)*, y_(i)*z_(i)*)∈Ω_(i). Itis assumed that the rate at which cells with attributes in Ω_(i) die attime t is proportional to the number of cells with attributes in Ω_(i)at time t, i.e., it is given by

α(t,x _(i) **,y _(i) **,z _(i)**)u(t,x _(i) *,y _(i) *,z _(i)*)

where (x_(i)**,y_(i)**,z_(i)**), (x_(i)*,y_(i)*,z_(i)*)ϵΩ_(i) anda(t,x,y,z) is the mortality rate, i.e., the rate at which cells withattributes (x,y,z) die at time t. Summing up and taking.epsilon..fwdarw.0, the rate at which individuals with attributes in Ωdie at time t is given by

∫_(Ω)α(t,x,y,z)u(t,x,y,z)dω.  (A.10)

c) Proliferation.

A basic assumption for proliferation is that division at time t of acell with cell age x, amount of hemoglobin y and amount of transferrinreceptors z results in two cells each with cell age x, amount ofhemoglobin y/2 and amount of transferrin receptors z/2. Furthermore, therate at which cells with attributes (x,y,z) divide at time t is assumedto be given by β(t,x,y,z).

For a subset Ω⊂R₊ ³ the subsets Ω^((1/2)), Ω⁽²⁾⊂R₊ ³ are defined by

${\Omega^{(\frac{1}{2})} = \left\{ {\left( {x,y,z} \right)❘{\left( {x,{2y},{2z}} \right) \in \Omega}} \right\}},$$\Omega^{2} = \left\{ {\left( {x,y,z} \right)❘{\left( {x,\frac{y}{2},\frac{z}{2}} \right) \in \Omega}} \right\}$

From now on, it is assumed that

Ω^((1/2))∩Ω=ϕ and Ω⁽²⁾∩Ω=ϕ.  (A.11)

The following results are important when the equations governing thedevelopment of a structured population are derived:

Lemma A.1 For any (x₀,y₀,z₀)ϵR₊ ³ with y₀>0 or z₀>0 there exists asubset Ω⊂R₊ ³ with (x₀, y₀, z₀)∈Ω satisfying condition (A.11).

Proof. The case y₀>0 is the only case considered. For

$\varepsilon = {\frac{y_{0}}{3} > 0}$

set

Ω={(x,y,z)∈R ₊ ³ ∥x−x ₀ |<ε,|y−y ₀ |<ε,|z−z ₀|<ε}

and assume that (x, y, z)∈Ω∩Ω^((1/2)). By the definition of Ω^((1/2)),this is equivalent to

|x−x ₀ |<ε,|y−y ₀|<ε,|2y−y ₀ |<ε,|z−z ₀|<ε, and |2z−z ₀|<ε

The inequalities for y and the definition of ε imply that y>y₀−ε=2y₀/3and y (y₀+ε)/2=2y₀/3, a contradiction, which proves that Ω^(1/2)∩Ω=Ø.The proof for Ω2∩Ω=Ø is analogous. The rate of change of the number ofcells with attributes in Ω due to proliferation equals the rate at whichthe daughter cells of cells with attributes in Ω₂ move into Ω the rateat which the daughter cells of cells with attributes in Ω leave Ω (andmove to Ω^((1/2))). By analogous considerations as in the case of theeffect of mortality, the following expression is obtained for the rateat which the number of cells with attributes in Ω. changes because ofproliferation:

$\begin{matrix}{{{2{\int_{\Omega^{(2)}}{{\beta\left( {t,\overset{\sim}{x},\overset{\sim}{y},\overset{\sim}{z}} \right)}d\overset{\sim}{\omega}}}} - {\int_{\Omega}{{\beta\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}d\omega}}} = {{8{\int_{\Omega}{{\beta\left( {t,x,{2y},{2z}} \right)}{u\left( {t,x,{2y},{2z}} \right)}d\omega}}} - {\int_{\Omega}{{\beta\left( {t,x,y,z} \right)}{u\left( {t,x,y,z} \right)}d{\omega.}}}}} & \left( {A.12} \right)\end{matrix}$

d) The Model.

Using (A.6), (A.9), (A.10) and (A.12) yields

∫_(Ω)u_(t)(t, x, y, z)dω = −∫_(Ω)((f(t, x, y, z)u(t, x, y, z))_(x) + (g(t, x, y, z)u(t, x, y, z)_(y) + (h(t, x, y, z)u(t, x, y, z))_(z))dω − ∫_(Ω)(α(t, x, y, z) + β(t, x, y, z))u(t, x, y, z)dω + 8∫_(Ω)β(t, x, 2y, 2z)u(t, x, 2y, 2z)dω.

Assuming that the integrands in the above integrals are continuous andusing Lemma 1 yields (A.13)

u₁(t,x,y,z)+(ƒ(t,x,y,z)u(t,x,y,z))_(x)+(g(t,x,y,z)u(t,x,y,z))_(y)+(h(t,x,y,z)u(t,x,y,z))_(z)=−(α(t,x,y,z)+β(t,x,y,z))u(t,x,y,z)+8β(t,x,2y,2z)u(t,x,2y,2z),t≥0,x>0,y>0,z>0.

An initial condition is prescribed

u(0,x,y,z)=u _(o)(x,y,z),x>0,y>0,z>0

and boundary conditions which, because of their special form, areconsidered directly in Section 2.A.3. List of Parameters for models described in Sections 6 and 7

TABLE 5 List of parameters and their meaning. Parameter Meaning β^(p)proliferation rate for BFU-E cells β^(q) proliferation rate for CFU-Ecells β^(r) proliferation rate for erythroblasts μ_(max) ^(p) maximalmaturity for BFU-E cells μ_(min) ^(q) minimal maturity for CFU-E cellsμ_(max) ^(q) maximal maturity for CFU-E cells μ_(min) ^(r) minimalmaturity for erythroblasts μ_(max) ^(r) maximal maturity forerythroblasts μ_(min)

minimal maturity for marrow reticulocytes μ_(max)

maximal maturity for marrow reticulocytes α

^(m) intrinsic mortality rate for erythrocytes a₁, b₁, c₁, k₁ constantsfor the sigmoid apoptosis rate for CFU-E cells a₂, b₂, c₂, k₂ constantsfor the sigmoid maturation velocity for marrow reticulocytes d₁ rateconstant for transferrin to bind to its receptor d rate constant withwhich hemoglobin is synthesized after transferrin bound to the TfR a₃,b₃, c₃, k₃ constants for the sigmoid function governing the release ofendogenous EPO μ_(min) ^(m, n) lower bound of erythrocytes which arepossibly exposed to neocytolysis μ_(max) ^(m, n) lower bound oferythrocytes which are possibly exposed to neocytolysis μ_(max) maximallife span for erythrocytes c_(E) constant in the mortality rate forerythrocytes τ_(E) EPO threshold for neocytolysis c_(deg) ^(end)degradation rate of EPO released by the kidneys c_(deg) ^(ex)degradation rate of administered EPO (Epoetin alfa) S_(o) number ofcells committing to the erythroid lineage TBV total blood volume a₄, b₄,c₄, k₄ constants for the polynomial function describing the number ofTfR expressed on erythroblasts η calendric age of the cell η^(s)calendric age of the reticulocyte c

constant in the function describing the apoptosis rate for cells withage η at time t h

reference hemoglobin value k

transfer rates from iron compartment i to iron compartment j a_(iv) theamount of iron administered intravenously a

the amount of iron in the duodenum a₅, b₅, c₅, k₅ constants for thesigmoid function governing the flow rate from the duodenum into theplasma compartment a_(H), b_(H), c_(H), k_(H) constants for the sigmoidfunction governing the release of hepcidin by the liver c

, c

, c_(i) constants for the auxiliary equation used in the sigmoidfunction governing the release of hepcidin by the liver c_(deg) ^(H)degradation rate of hepcidin a₆, b₆, c₆, k₆ constants for the sigmoidfunction describing the flow rate from the macrophages into the plasmacompartment a₇, b₇, c₇, k₇ constants for the sigmoid function governingthe flow rate from the storage into the plasma compartment

indicates data missing or illegible when filed

A.4. Derivation of the Hgb Accumulation Function.

Let h(t,μ) respectively Ø(p) denote the amount of hemoglobin in a cellof maturity μ≥μ_(min) at time tϵR respectively the number of transferrinreceptors on the membrane of a cell of maturity μ≥μ_(min). Then the rateat which h(t,μ) changes in time is given by

(A.4.1)

${{\frac{d}{dt}{h\left( {t,{\mu(t)}} \right)}} = {{{\frac{\partial}{\partial t}{h\left( {t,{\mu(t)}} \right)}} + {{\frac{\partial}{\partial t}{h\left( {t,{u(t)}} \right)}}{v(t)}}} = {d{\varnothing\left( {\mu(t)} \right)}{\alpha_{1}(t)}}}},{t \in R}$

where d>0 is a constant, α₁(t) denotes the iron concentration in plasmaat time t and v(t) is the velocity at which the maturity of a cellincreases, v(t)=μ(t). For equation (A.4.1) one has the boundarycondition

(A.4.2)

μ(t,μ _(min))≡0,t∈R

We solve the problem (A.4.1), (A.4.2) using the method ofcharacteristics. Let s→(t(s),μ(s),u(s)), x≥0 be the characteristicthrough a point (t₀,μ). Then the functions t(s), μ(s) and u(s) aresolutions of the following initial value problem (“prime(′)”=d/ds):

(A.4.3)

t′(s)=1,t(0)=t ₀,

u′(s)=v(t(s)),μ(0)=μ_(min),

u′(s)=dØ(μ(s))a ₁(t(s)),u(0)=0

This gives, for

s≥0,

t(s)=t ₀ +s,  (A.4.4)

μ(s)=μ_(min)+∫₀ ^(s) v(t ₀+σ)dσ,  (A.4.5)

u(s)=h(t(s),μ(s))=d∫ ₀ ^(s)Ø(μ_(min)+∫₀ ^(σ) v(t ₀+ρ)dρ)a ₁(t₀+σ)dσ  (A.4.6)

In order to determine h(t,μ) with t∈R and μ≥μ_(min) one has to computetwo numbers s≥0 and {circumflex over (t)}∈R such that t=t(s),μ=μ(s) andh(t,μ)=u(s) where t(⋅), .μ(⋅) and u(⋅) solve system (A.4.3) witht₀={circumflex over (t)} From equations (A.4.4) and (A.4.5) for s=s weget

${\overset{\hat{}}{t} = {\overset{¯}{t} - \overset{¯}{s}}},{\overset{¯}{\mu} = {{\mu_{\min} + {\int_{0}^{\overset{¯}{t} - \overset{\hat{}}{t}}{{\nu\left( {\overset{\hat{}}{t} + \sigma} \right)}d\sigma}}} = {\mu_{m\overset{˙}{i}n} + {\int_{\hat{t}}^{\overset{¯}{t}}{{v(\sigma)}d{\sigma.}}}}}}$

This implies s=t−t and that {circumflex over (t)} has to satisfy theequation

(A.4.7)

${g\left( \overset{\hat{}}{t} \right)} = {{\int_{\overset{\hat{}}{t}}^{\overset{¯}{t}}{{v(\sigma)}d\sigma}} = {\overset{\_}{\mu} - \mu_{\min}}}$

The function {circumflex over (t)}→g({circumflex over (t)}) ismonotonically decreasing with g(t)=0. Equation (A.4.7) has a uniquesolution {circumflex over (t)}<t if, for instance, v(t)≥δ for t∈[T,{circumflex over (t)}] with δT≥μ_(max)−μ_(min) or if for some t*∈R onehas

∫_(−∞)^(t^(*))v(σ)dσ = ∞

In case

v≡l one gets μ=μ _(min)+t−{circumflex over (t)} and s=t−{circumflex over(t)}, i.e., {circumflex over (t)}=t−(μ−μ_(min)) and s=μ−μ_(min).

Furthermore, one has

${h\left( {\overset{¯}{t},\overset{¯}{\mu}} \right)} = {{u\left( \overset{¯}{s} \right)} = {{d{\int_{0}^{\overset{\_}{\mu} - \mu_{\min}}{{\varnothing\left( {\mu_{\min} + \sigma} \right)}{a_{1}\left( {\sigma + \overset{\_}{t} - \left( {\overset{¯}{\mu} - \mu_{\min}} \right)} \right)}d\sigma}}} = {d{\int_{- {({\overset{\_}{\mu} - \mu_{\min}})}}^{0}{{\varnothing\left( {\overset{\_}{\mu} + \sigma} \right)}{a_{1}\left( {\overset{\_}{t} + \sigma} \right)}d\sigma}}}}}$

Anemia of Renal Disease

Anemia affects almost all patients with chronic kidney disease (CKD). Itis caused by failure of renal excretory and endocrine function. It oftendevelops once renal function decreases to 50% and the degree of anemiaincreases with severity of renal failure. See Strippoli, et. al. Anemiadevelops because of a deficiency in endogenous erythropoietin productionby the kidneys, increased blood losses, functional or absolute irondeficiency and decreased red cell survival. See J. B. Wish, ClinicalJournal of the American Society of Nephrology, 1:S4-S8, 2006. Concerningendogenous erythropoietin production by the kidneys, one has to notethat kidneys of end-stage renal disease patients can still produce someEPO, and thus can maintain higher hemoglobin levels than those found inanephric patients. Reasons for blood losses can be purpura,gastrointestinal and gynecologic bleeding (which occur in one third toone half of all CKD patients), frequent blood sampling, blood left inthe extracorporeal circuit, multiple vascular surgeries, etc.Furthermore, neocytolysis contributes to renal anemia and it explainsthe often demonstrable hemolytic component and the worsening ofhemolysis with more pronounced renal disease. A description ofneocytolysis is provided above. It further explains the responsivenessof hemolysis to ESA therapy and the increased efficiency of subcutaneous(s.c.) compared to intravenous therapy (i.v.), because in s.c.administration nadirs in EPO levels which precipitate neocytolysis areavoided. See Rice 1999. In general, renal anemia is normocytic andnormochromic and the number of reticulocytes is normal or slightlydecreased, which is inappropriate in the context of a reduced RBCpopulation. Certain deficiency states, especially iron, but also folateor vitamin B.sub.12 deficiency may alter the nature of the anemia. SeeFinch 1982; M. Polenakovic and A. Sikole, Journal of American Society ofNephrology, 7:1178-1182, 1996.

Untreated anemia is associated with decreased oxygen delivery to thetissues. For compensatory reasons, cardiac output increases, resultingin left ventricular hypertrophy. Cardiac disease is the most commoncause of death among patients who are on maintenance dialysis. Partialcorrection of anemia in these patients was shown to reduce cardiacischemia and ameliorate cardiomegaly, thus reducing cardiac relatedmorbidity. See A. Besarab, W. K. Bolton, J. K. Browne, J. C. Egrie, A.R. Nissenson, D. M. Okamoto, S. J. Schwab, and D. A. Goodkin, The NewEngland Journal of Medicine, 339:584-590, 1998. Further consequences ofuncorrected anemia are decreased cognition and mental acuity and overalldecrease in patient welfare.

ESA Therapy

Erythropoiesis stimulating agents have been used to treat anemia inpatients suffering from chronic renal failure for more than two decades.Optimal hemoglobin targets are still a matter of discussion. Studieshave shown (see e.g. Singh et al.; Strippoli et al., R. N. Foley, B. M.Curtis, and P. S. Parfrey. Clinical Journal of the American Society ofNephrology, 4:726-733, 2009 (hereinafter “Foley et al.”) that a partialcorrection of anemia is preferable to a full correction. A number ofcomplications can occur in patients with CKD when they havenear-normal/normal hemoglobin levels, i.e., higher vascular accessthrombosis, hypertension and greater requirements for antihypertensives,cardiovascular events, earlier need for renal replacement therapy andhigher mortality. Normal hemoglobin levels for women are in a range ofbetween about 12 g/dl and 16 g/dl, and in a range of between about 13g/dl and about 17.5 g/dl for men.

Defining an optimal hemoglobin target is not the only issue regardingESA therapy. Another problem is: how can the target hemoglobin beachieved and how to keep the patient near this hemoglobin level over along time period? The dose and frequency of administration of an ESAtreatment regimen are most often determined based on the priorexperience of the physician and on established guidelines. This approachbears some limitations and level of hemoglobin tends to fluctuategreatly and cycling phenomena are observed. An analysis of 31,267patients on hemodialysis in the Fresenius Medical Care-North Americadatabase found that only 5% of patients persistently remained within adesired Hb range of 11-12 g/dl for a period of 6 months. See A. J.Collins, R. M. Brenner, J. J. Ofman, E. M. Chi, N. Stuccio-White, M.Krishnan, C. Solid, N. J. Ofsthun, and J. M. Lazarus. American Journalof Kidney Diseases, 46:481-488, 2005 (hereinafter “Collins et al.”).Fishbane et al. analyzed data of dialysis patients collected over fiveyears in a hospital and came to a similar conclusion. See S. Fishbaneand J. S. Berns. Kidney International, 68:1337-1343, 2005 (hereinafter“Fishbane et al.”). More than 90% of the patients experienced hemoglobincycling. The authors state that changes in ESA dose were the mostimportant driver and were associated with hemoglobin excursion in about80% of cases. The ESA dose adjustment protocol that was used was similarto the protocol used in most dialysis centers.

Therapy with ESAs is quite different than biological erythropoietinsecretion. In hemodialysis patients, i.v. administration route isprimarily used, because of the availability of venous access.Intravenous treatment involves short, intermittent, non-physiologicbursts of EPO concentration. The bursts are followed by a fast declineto very low levels of EPO. These fluctuations in plasma concentration donot coincide, either temporally or in magnitude, with physiologicperturbations. Therefore, it may not be surprising that Hb levelsfluctuate widely and that it is extremely difficult for physicians toadjust dosing schemes such that no cycling phenomena occur. Note thatfluctuations in hemoglobin result in varying oxygen delivery to vitalorgans. Consequences include repeated episodes of relative ischemia andcompensatory mechanisms in organs (e.g., heart) that may result indisordered growth signals, pathologic organ function and worsenedpatient outcomes. (See Fishbane et al.).

A predictive model of erythropoiesis can help deal with this situation.For instance, a physician can try different ESA treatment regimens andobserve their effects on hemoglobin levels over the next few months.Dosing regimens can be tested/chosen with regard to avoidance ofneocytolysis, minimal amounts of ESA administered and avoidance ofcycling patterns in Hb concentration.

Iron Therapy

Patients with anemia of chronic kidney disease have to be followed forsymptoms of iron deficiency. See J. W. Fisher. Experimental Biology andMedicine, 28:1-24, 2003. 80-90% of dialysis patients on ESA therapy willrequire iron at some stage. See R. M. Schaefer and L. Schaefer.Nephrology Dialysis Transplantation, 13:9-12, 1998. This very pronouncedneed for supplementary iron has different reasons. On one hand ironstores can be depleted (absolute iron deficiency). Iron loss inhemodialysis patients (due to continuous gastrointestinal, purpural andgynecological bleeding, frequent blood sampling, surgeries, . . . ) isabout 1500-3000 mg/year, as compared to iron loss for a healthy adultthat is about 400-800 mg/year, and can be even higher under certaincircumstances. Hence, the daily need for iron can be well above theabsorptive capacity of the intestinal. This is aggravated by the factthat the uptake via the duodenum is often impaired in these patients. Onthe other hand, a functional iron deficiency is often observed in renalanemia. During ESA therapy, the number of erythroid progenitor andprecursor cells in the bone marrow increase drastically and this imposesa lot of stress on systemic iron homeostasis. Supply and demand often donot match. It is difficult, even for healthy persons, to increase therate at which iron is released from stores and is recycled fromhemoglobin to deliver enough iron to the bone marrow to keep up withsupraphysiologic rates of RBC production during ESA treatment. Mattersare further complicated in chronic kidney disease because ofinflammation, which frequently occurs with various degrees of severity.Thus, iron utilization is regularly decreased in renal insufficiency.See above for a detailed description of the effects of inflammation onsystemic iron homeostasis.

In the pre-ESA era iron overload, because of excessive bloodtransfusions, was a major cause of morbidity in dialysis patients. Itssignificance in the post-EPO era remains unclear. See K. Kalantar-Zadeh,D. L. Regidor, C. J. McAllister, B. Michael, and Warnock D. G., Journalof American Society of Nephrology, 16:3070-3080, 2005. A continuousadministration of i.v. iron certainly results in overfilled iron-storesand imposes serious health risks. Hence, decisions on when to supply thepatient with i.v. iron and when to withdraw therapy have to bethoroughly evaluated. A mathematical model can help keep track of thecurrent iron status of a patient and can help to make decisions onadaptations of treatment.

Finally, there is a very complex interaction between hemoglobin cyclingand iron storage whose dynamical behavior under ESA and iron treatment,even for an experienced physician, is barely predictable. In Alfrey etal., the authors suggest that: “[ . . . ] the current therapeuticparadigm of hemoglobin monitoring, iron treatment, and rHuEPO treatmentresults in recurrent nonphysiologic cycling of hemoglobin levels inhemodialysis patients.”

Computation of Hematocrit and Hemoglobin Concentrations

In order to compute the hematocrit (HCT) and the hemoglobinconcentration (Hb) for a subject from the model output (which is thenumber of red blood cells circulating in blood), estimates of the totalblood volume of the subject are needed. The Nadler and Allen formula(see Nadler, S. B., Hidalgo, J. U., and Bloch, T., Surgery, 1962, Vol.51, pp. 224-232) can be used to estimate the total blood volume (TBV) ofa healthy subject according to her/his weight and height:

TBV [ml]=183.3+356.1·times·(height [m])³+33.08·times·weight[kg],  Women:

TBV [ml]=604.1+366.9·times·(height [m])³+32.19·times·weight [kg].  Men:

The TBV for a dialysis patient can be measured by radio-labeling redblood cells with chromium-51 as described above. Using these estimatesfor the total blood volume, the hematocrit (HCT) and the hemoglobinconcentration (Hb) can be computed for a patient from the number of redblood cells circulating in blood via the following formulae.

${HC{T\lbrack\%\rbrack}} = \frac{\left( {{M(t)} \times {{MCV}\lbrack{fl}\rbrack}} \right)}{{TBV}\lbrack{ml}\rbrack}$

where M(t)=∫₀ ^(μ) ^(max) ^(m) m(t, μ^(m))dμ^(m) is the total number oferythrocytes circulating in blood and MCV is the mean corpuscular volumeof a RBC which is obtained from measurement, a

${H{b\left\lbrack \frac{g}{l} \right\rbrack}} = {1000 \times \frac{\left( {{M(t)} \times {{MCV}\left\lbrack {pg} \right\rbrack}} \right)}{{TBV}\lbrack{ml}\rbrack}}$

where MCH is the mean cellular hemoglobin, which is also obtained viameasurements.

Implementation in a Computer Network

FIG. 14 illustrates a computer network or similar digital processingenvironment in which the present invention can be implemented.

Computer(s)/devices 50 and server computer(s) 60 provide processing,storage, and input/output devices executing application programs and thelike. Computer(s)/devices 50 can also be linked through communicationsnetwork 70 to other computing devices, including other devices/processes50, digital processor dialysis machines 50A, including machines withintegrated modular drug delivery devices, and server computer(s) 60.Communications network 70 can be part of a remote access network, aglobal network (e.g., the Internet), a worldwide collection ofcomputers, Local area or Wide area networks, and gateways that currentlyuse respective protocols (TCP/IP, Bluetooth, etc.) to communicate withone another. Other electronic device/computer network architectures aresuitable.

FIG. 15 is a diagram of the internal structure of a computer (e.g.,processor/device 50, digital processor dialysis machines 50A, or servercomputers 60) in the computer system of FIG. 14 . Each computer 50, 50A,60 contains a system bus 79, where a bus is a set of hardware lines usedfor data transfer among the components of a computer or processingsystem. Bus 79 is essentially a shared conduit that connects differentelements of a computer system (e.g., processor, disk storage, memory,input/output ports, network ports, etc.) that enables the transfer ofinformation between the elements which could be data collected during atreatment, but also instructions to a machine with integrated modulardrug delivery devices to follow a drug delivery scheme provided by acentral erythropoiesis modeling engine. Attached to system bus 79 is I/Odevice interface 82 for connecting various input and output devices(e.g., keyboard, mouse, displays, printers, speakers, etc.) to thecomputer 50, 60. Network interface 86 allows the computer to connect tovarious other devices attached to a network (e.g., network 70 of FIG. 14). Memory 90 provides volatile storage for computer softwareinstructions 92 and data 94 used to implement an embodiment of thepresent invention (e.g., Eqs. described above and listed in Section 4 orany other erythropoiesis modeling engine detailed above). Disk storage95 provides non-volatile storage for computer software instructions 92and data 94 used to implement an embodiment of the present invention.Central processor unit 84 is also attached to system bus 79 and providesfor the execution of computer instructions.

In one embodiment, the processor routines 92 and data 94 are a computerprogram product (generally referenced 92), including a computer readablemedium (e.g., a removable storage medium such as one or more DVD-ROM's,CD-ROM's, diskettes, tapes, etc.) that provides at least a portion ofthe software instructions for the system. Computer program product 92can be installed by any suitable software installation procedure, as iswell known in the art. In another embodiment, at least a portion of thesoftware instructions can also be downloaded over a cable, communicationand/or wireless connection. In other embodiments, the programs are acomputer program propagated signal product 107 embodied on a propagatedsignal on a propagation medium (e.g., a radio wave, an infrared wave, alaser wave, a sound wave, or an electrical wave propagated over a globalnetwork such as the Internet, or other network(s)). Such carrier mediumor signals provide at least a portion of the software instructions forthe present routines/program 92.

In alternate embodiments, the propagated signal is an analog carrierwave or digital signal carried on the propagated medium. For example,the propagated signal may be a digitized signal propagated over a globalnetwork (e.g., the Internet), a telecommunications network, or othernetwork. In one embodiment, the propagated signal is a signal that istransmitted over the propagation medium over a period of time, such asthe instructions for a software application sent in packets over anetwork over a period of milliseconds, seconds, minutes, or longer. Inanother embodiment, the computer readable medium of computer programproduct 92 is a propagation medium that the computer system 50 canreceive and read, such as by receiving the propagation medium andidentifying a propagated signal embodied in the propagation medium, asdescribed above for computer program propagated signal product.

Generally speaking, the term “carrier medium” or transient carrierencompasses the foregoing transient signals, propagated signals,propagated medium, storage medium and the like.

The relevant teachings of all patents, published patent applications andliterature references cited herein are incorporated by reference intheir entirety.

While this invention has been particularly shown and described withreferences to example embodiments thereof, it will be understood bythose skilled in the art that various changes in form and details may bemade therein without departing from the scope of the inventionencompassed by the appended claims.

The meanings of method steps of the invention(s) described herein areintended to include any suitable method or manner of causing one or moreother parties or entities to perform the recited steps, consistent withthe patentability of the following claims, unless a different meaning isexpressly provided or otherwise clear from the context. Such parties orentities need not be under the direction or control of any other partyor entity, and need not be located within a particular jurisdiction.

Thus for example, a description or recitation of “adding a first numberto a second number” is meant to include causing one or more parties orentities to add the two numbers together. For example, if person Xengages in an arm's length transaction with person Y to add the twonumbers, and person Y indeed adds the two numbers, then both persons Xand Y perform the step as recited: person Y by virtue of the fact thathe actually added the numbers, and person X by virtue of the fact thathe caused person Y to add the numbers. Furthermore, if person X islocated within the United States and person Y is located outside theUnited States, then the method is performed in the United States byvirtue of person X's participation in causing the step to be performed.

1. A method for determining efficacy of a therapy, comprising: a)generating in silico a plurality of virtual patients for modeling ahealth condition based on data collected from a population of previouslytreated patients; wherein the collected data represents at least onemeasured biological response of said previously treated patients to apreviously administered therapeutic regimen; b) applying a simulatedtherapy to said virtual patients; and c) determining one or morephysiological parameters in said virtual patients in response to saidsimulated therapy.
 2. The method of claim 1, further comprising: d)adjusting said simulated therapy based on said determined one or morephysiological parameters of said virtual patients; e) applying saidadjusted simulated therapy to said virtual patients; and f) determiningthe one or more physiological parameters in response to said adjustedsimulated therapy.
 3. The method of claim 2, further comprisingrepeating steps d)-f) so as to optimize said simulated therapy forapplication to one or more actual patients.
 4. The method of claim 1,further comprising determining whether at least one of said determinedphysiological parameters of at least one of said virtual patients isindicative of an adverse effect from said simulated therapy.
 5. Themethod of claim 1, wherein the simulated therapy includes a simulatednoncompliance with a prescribed therapy.
 6. The method of claim 1,wherein the determining step further comprises determining values ofsaid one or more physiological parameters at a plurality of times over apredetermined time interval.
 7. The method of claim 6, wherein the oneor more physiological parameters comprises one or more metabolicparameters.
 8. The method of claim 7, wherein the one or more metabolicparameters comprise blood parameters and/or urine parameters.
 9. Themethod of claim 8, wherein the blood parameters comprise any ofhemoglobin concentration and hematocrit level.
 10. The method of claim1, wherein each of the virtual patients comprises at least one in silicomodel representing a physiological system.
 11. The method of claim 10,wherein the at least one model is an erythropoiesis model.
 12. Themethod of claim 1, further comprising applying a plurality of differentsimulated therapies to each of said virtual patients, determining one ormore physiological parameters in said virtual patients in response toeach simulated therapy, and evaluating each simulated therapy based onthe corresponding physiological parameters.
 13. The method of claim 1,wherein the therapy comprises any of at least one of: a pharmacologicaltherapy, a non-pharmacological therapy, and a combination thereof 14.The method of claim 13, wherein the at least one pharmacological therapycomprises an anemia therapy.
 15. The method of claim 14, wherein theanemia therapy comprises administering at least one of an erythropoiesisstimulating agent (ESA), iron, a drug stimulating endogenouserythropoietin release and/or synthesis, a biosimilar, or a combinationthereof.
 16. The method of claim 15, wherein the erythropoiesisstimulating agent comprises exogenous erythropoietin.
 17. The method ofclaim 13, wherein the one or more non-pharmacological therapiescomprises a fluid therapy, a dietary therapy, an exercise therapy, anextracorporeal therapy, a radiotherapy, a therapy using sound and/orultrasound, an electrotherapy, or a combination thereof
 18. The methodof claim 1, wherein the data collected from the population of previouslytreated patients comprises gender, age, weight, height, ethnicity,metabolic/chemistry parameters, complete blood count, or a combinationthereof.
 19. The method of claim 18, wherein the data collected furthercomprises medications used by the previously treated patients, pastmedical history data, past surgical history data or a combinationthereof
 20. The method of claim 19, wherein the past medical historydata comprises information regarding diabetes, bloodpressure/hypertension, cancer, congestive heart failure, or acombination thereof. 21-54. (canceled)